Stokes I Simultaneous Image and Instrument Modeling

In this tutorial, we will create a preliminary reconstruction of the 2017 M87 data on April 6 by simultaneously creating an image and model for the instrument. By instrument model, we mean something akin to self-calibration in traditional VLBI imaging terminology. However, unlike traditional self-cal, we will at each point in our parameter space effectively explore the possible self-cal solutions. This will allow us to constrain and marginalize over the instrument effects, such as time variable gains.

To get started we load Comrade.

using Comrade
  Activating project at `~/work/Comrade.jl/Comrade.jl/examples`
using Pyehtim
using LinearAlgebra

For reproducibility we use a stable random number genreator

using StableRNGs
rng = StableRNG(123)
StableRNGs.LehmerRNG(state=0x000000000000000000000000000000f7)

Load the Data

To download the data visit https://doi.org/10.25739/g85n-f134 First we will load our data:

obs = ehtim.obsdata.load_uvfits(joinpath(dirname(pathof(Comrade)), "..", "examples", "SR1_M87_2017_096_hi_hops_netcal_StokesI.uvfits"))
Python: <ehtim.obsdata.Obsdata object at 0x7f6d17113af0>

Now we do some minor preprocessing:

  • Scan average the data since the data have been preprocessed so that the gain phases coherent.
  • Add 1% systematic noise to deal with calibration issues that cause 1% non-closing errors.
obs = scan_average(obs.add_fractional_noise(0.02))
Python: <ehtim.obsdata.Obsdata object at 0x7f6d366b4d00>

Now we extract our complex visibilities.

dvis = extract_table(obs, ComplexVisibilities())
EHTObservation{Float64,Comrade.EHTVisibilityDatum{Float64}, ...}
  source: M87
  mjd: 57849
  frequency: 2.29070703125e11
  bandwidth: 1.856e9
  stations: [:AA, :AP, :AZ, :JC, :LM, :PV, :SM]
  nsamples: 284

##Building the Model/Posterior

Now, we must build our intensity/visibility model. That is, the model that takes in a named tuple of parameters and perhaps some metadata required to construct the model. For our model, we will use a raster or ContinuousImage for our image model. The model is given below:

function sky(θ, metadata)
    (;fg, c, σimg) = θ
    (;ftot, K, meanpr, grid, cache) = metadata
    # Transform to the log-ratio pixel fluxes
    cp = meanpr .+ σimg.*c.params
    # Transform to image space
    rast = (ftot*(1-fg))*K(to_simplex(CenteredLR(), cp))
    img = IntensityMap(rast, grid)
    m = ContinuousImage(img, cache)
    # Add a large-scale gaussian to deal with the over-resolved mas flux
    g = modify(Gaussian(), Stretch(μas2rad(250.0), μas2rad(250.0)), Renormalize(ftot*fg))
    return m + g
end
sky (generic function with 1 method)

Unlike other imaging examples (e.g., Imaging a Black Hole using only Closure Quantities) we also need to include a model for the instrument, i.e., gains as well. The gains will be broken into two components

  • Gain amplitudes which are typically known to 10-20%, except for LMT, which has amplitudes closer to 50-100%.
  • Gain phases which are more difficult to constrain and can shift rapidly.
function instrument(θ, metadata)
    (; lgamp, gphase) = θ
    (; gcache, gcachep) = metadata
    # Now form our instrument model
    gvis = exp.(lgamp)
    gphase = exp.(1im.*gphase)
    jgamp = jonesStokes(gvis, gcache)
    jgphase = jonesStokes(gphase, gcachep)
    return JonesModel(jgamp*jgphase)
end
instrument (generic function with 1 method)

The model construction is very similar to Imaging a Black Hole using only Closure Quantities, except we include a large scale gaussian since we want to model the zero baselines. For more information about the image model please read the closure-only example. Let's discuss the instrument model Comrade.JonesModel. Thanks to the EHT pre-calibration, the gains are stable over scans. Therefore, we can model the gains on a scan-by-scan basis. To form the instrument model, we need our

  1. Our (log) gain amplitudes and phases are given below by lgamp and gphase
  2. Our function or cache that maps the gains from a list to the stations they impact gcache.
  3. The set of Comrade.JonesPairs produced by jonesStokes

These three ingredients then specify our instrument model. The instrument model can then be combined with our image model cimg to form the total JonesModel.

Now, let's set up our image model. The EHT's nominal resolution is 20-25 μas. Additionally, the EHT is not very sensitive to a larger field of view. Typically 60-80 μas is enough to describe the compact flux of M87. Given this, we only need to use a small number of pixels to describe our image.

npix = 32
fovx = μas2rad(150.0)
fovy = μas2rad(150.0)
7.27220521664304e-10

Now let's form our cache's. First, we have our usual image cache which is needed to numerically compute the visibilities.

grid = imagepixels(fovx, fovy, npix, npix)
buffer = IntensityMap(zeros(npix, npix), grid)
cache = create_cache(NFFTAlg(dvis), buffer, BSplinePulse{3}())
VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}(VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}(NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}(1, 4, 2.0, :kaiser_bessel, AbstractNFFTs.POLYNOMIAL, true, false, true, 0x00000000), [1.2065151329523807e9 819926.8691406249 … 4.472903950222221e9 -15688.352240668402; -3.8185738483809524e9 -1.7123166484375e6 … -6.2063795222222224e7 119336.55902777778]), NFFTPlan with 284 sampling points for an input array of size(32, 32) and an output array of size(284,) with dims 1:2, ComplexF64[0.9304270810105981 - 0.1755525366839422im, 0.9999999857224657 - 6.37119066797087e-5im, 0.6968566352375541 - 0.5221791050707502im, 0.9304860627318433 + 0.17550202040945986im, 0.8028821861479044 + 0.3950138433443315im, 0.696909133455242 + 0.522149977275377im, 0.9999999859666728 - 6.454986502401626e-5im, 0.690100908057055 - 0.5259755617707915im, 0.9317353097306306 - 0.1646343527093735im, 0.9317937328057632 + 0.1645830409592267im  …  0.8346090489457091 + 0.1795789988845693im, 0.834610800886739 + 0.17958571977595708im, 0.9181402138641586 + 0.09206131734480663im, 0.9120864418553256 + 0.2942607208902332im, 0.9120838428437319 + 0.29426742718619836im, 0.9559925804493471 + 0.09847981391058878im, 0.8881131676530709 + 0.289296040490046im, 0.9924358388500316 + 0.0028007013218525043im, 0.8881114735852544 + 0.28930284794286826im, 0.9999999999233905 + 7.3999334332420275e-6im], BSplinePulse{3}(), [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0])

Second, we now construct our instrument model cache. This tells us how to map from the gains to the model visibilities. However, to construct this map, we also need to specify the observation segmentation over which we expect the gains to change. This is specified in the second argument to jonescache, and currently, there are two options

  • FixedSeg(val): Fixes the corruption to the value val for all time. This is usefule for reference stations
  • ScanSeg(): which forces the corruptions to only change from scan-to-scan
  • TrackSeg(): which forces the corruptions to be constant over a night's observation

For this work, we use the scan segmentation for the gain amplitudes since that is roughly the timescale we expect them to vary. For the phases we use a station specific scheme where we set AA to be fixed to unit gain because it will function as a reference station.

gcache = jonescache(dvis, ScanSeg())
gcachep = jonescache(dvis, ScanSeg(); autoref=SEFDReference((complex(1.0))))

using VLBIImagePriors

Now we need to specify our image prior. For this work we will use a Gaussian Markov Random field prior Since we are using a Gaussian Markov random field prior we need to first specify our mean image. This behaves somewhat similary to a entropy regularizer in that it will start with an initial guess for the image structure. For this tutorial we will use a a symmetric Gaussian with a FWHM of 60 μas

fwhmfac = 2*sqrt(2*log(2))
mpr = modify(Gaussian(), Stretch(μas2rad(50.0)./fwhmfac))
imgpr = intensitymap(mpr, grid)
2-dimensional KeyedArray(NamedDimsArray(...)) with keys:
↓   X ∈ 32-element LinRange{Float64,...}
→   Y ∈ 32-element LinRange{Float64,...}
And data, 32×32 NamedDimsArray(::Matrix{Float64}, (:X, :Y)):
                 (-3.52247e-10)  …  (3.29522e-10)  (3.52247e-10)
 (-3.52247e-10)     6.37537e-8        1.32434e-7     6.37537e-8
 (-3.29522e-10)     1.32434e-7        2.751e-7       1.32434e-7
 (-3.06796e-10)     2.62014e-7        5.44273e-7     2.62014e-7
 (-2.84071e-10)     4.93724e-7        1.0256e-6      4.93724e-7
 (-2.61345e-10)     8.86092e-7   …    1.84065e-6     8.86092e-7
 (-2.38619e-10)     1.51463e-6        3.14629e-6     1.51463e-6
    ⋮                            ⋱    ⋮            
  (2.15894e-10)     2.46586e-6        5.12225e-6     2.46586e-6
  (2.38619e-10)     1.51463e-6        3.14629e-6     1.51463e-6
  (2.61345e-10)     8.86092e-7   …    1.84065e-6     8.86092e-7
  (2.84071e-10)     4.93724e-7        1.0256e-6      4.93724e-7
  (3.06796e-10)     2.62014e-7        5.44273e-7     2.62014e-7
  (3.29522e-10)     1.32434e-7        2.751e-7       1.32434e-7
  (3.52247e-10)     6.37537e-8        1.32434e-7     6.37537e-8

Now since we are actually modeling our image on the simplex we need to ensure that our mean image has unit flux

imgpr ./= flux(imgpr)
2-dimensional KeyedArray(NamedDimsArray(...)) with keys:
↓   X ∈ 32-element LinRange{Float64,...}
→   Y ∈ 32-element LinRange{Float64,...}
And data, 32×32 Matrix{Float64}:
                 (-3.52247e-10)  …  (3.29522e-10)  (3.52247e-10)
 (-3.52247e-10)     6.38049e-8        1.3254e-7      6.38049e-8
 (-3.29522e-10)     1.3254e-7         2.75321e-7     1.3254e-7
 (-3.06796e-10)     2.62224e-7        5.4471e-7      2.62224e-7
 (-2.84071e-10)     4.94121e-7        1.02642e-6     4.94121e-7
 (-2.61345e-10)     8.86803e-7   …    1.84213e-6     8.86803e-7
 (-2.38619e-10)     1.51585e-6        3.14882e-6     1.51585e-6
    ⋮                            ⋱    ⋮            
  (2.15894e-10)     2.46784e-6        5.12636e-6     2.46784e-6
  (2.38619e-10)     1.51585e-6        3.14882e-6     1.51585e-6
  (2.61345e-10)     8.86803e-7   …    1.84213e-6     8.86803e-7
  (2.84071e-10)     4.94121e-7        1.02642e-6     4.94121e-7
  (3.06796e-10)     2.62224e-7        5.4471e-7      2.62224e-7
  (3.29522e-10)     1.3254e-7         2.75321e-7     1.3254e-7
  (3.52247e-10)     6.38049e-8        1.3254e-7      6.38049e-8

and since our prior is not on the simplex we need to convert it to unconstrained or real space.

meanpr = to_real(CenteredLR(), Comrade.baseimage(imgpr))
32×32 Matrix{Float64}:
 -7.55422  -6.82317  -6.14085  -5.50727   …  -6.14085  -6.82317  -7.55422
 -6.82317  -6.09211  -5.4098   -4.77622      -5.4098   -6.09211  -6.82317
 -6.14085  -5.4098   -4.72748  -4.0939       -4.72748  -5.4098   -6.14085
 -5.50727  -4.77622  -4.0939   -3.46032      -4.0939   -4.77622  -5.50727
 -4.92243  -4.19137  -3.50906  -2.87548      -3.50906  -4.19137  -4.92243
 -4.38632  -3.65527  -2.97295  -2.33937   …  -2.97295  -3.65527  -4.38632
 -3.89895  -3.1679   -2.48558  -1.852        -2.48558  -3.1679   -3.89895
 -3.46032  -2.72927  -2.04695  -1.41337      -2.04695  -2.72927  -3.46032
 -3.07043  -2.33937  -1.65705  -1.02348      -1.65705  -2.33937  -3.07043
 -2.72927  -1.99821  -1.3159   -0.682317     -1.3159   -1.99821  -2.72927
  ⋮                                       ⋱             ⋮        
 -3.07043  -2.33937  -1.65705  -1.02348      -1.65705  -2.33937  -3.07043
 -3.46032  -2.72927  -2.04695  -1.41337      -2.04695  -2.72927  -3.46032
 -3.89895  -3.1679   -2.48558  -1.852     …  -2.48558  -3.1679   -3.89895
 -4.38632  -3.65527  -2.97295  -2.33937      -2.97295  -3.65527  -4.38632
 -4.92243  -4.19137  -3.50906  -2.87548      -3.50906  -4.19137  -4.92243
 -5.50727  -4.77622  -4.0939   -3.46032      -4.0939   -4.77622  -5.50727
 -6.14085  -5.4098   -4.72748  -4.0939       -4.72748  -5.4098   -6.14085
 -6.82317  -6.09211  -5.4098   -4.77622   …  -5.4098   -6.09211  -6.82317
 -7.55422  -6.82317  -6.14085  -5.50727      -6.14085  -6.82317  -7.55422

Now we can form our metadata we need to fully define our model.

metadata = (;ftot=1.1, K=CenterImage(imgpr), meanpr, grid, cache, gcache, gcachep)
(ftot = 1.1, K = VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}([0.9944957386363635 -0.005326704545454545 … 0.00532670454545453 0.0055042613636363605; -0.005326704545454545 0.9948393969941347 … 0.0051606030058651 0.005326704545454565; … ; 0.00532670454545453 0.0051606030058651 … 0.9948393969941348 -0.005326704545454544; 0.0055042613636363605 0.005326704545454565 … -0.005326704545454544 0.9944957386363638], (32, 32)), meanpr = [-7.55422122563378 -6.823167558636962 … -6.823167558636962 -7.55422122563378; -6.823167558636962 -6.092113891640146 … -6.092113891640146 -6.823167558636962; … ; -6.823167558636962 -6.092113891640146 … -6.092113891640146 -6.823167558636962; -7.55422122563378 -6.823167558636962 … -6.823167558636962 -7.55422122563378], grid = GriddedKeys{(:X, :Y)}
	X: LinRange{Float64}(-3.5224744018114725e-10, 3.5224744018114725e-10, 32)
	Y: LinRange{Float64}(-3.5224744018114725e-10, 3.5224744018114725e-10, 32)
, cache = VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}(VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}(NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}(1, 4, 2.0, :kaiser_bessel, AbstractNFFTs.POLYNOMIAL, true, false, true, 0x00000000), [1.2065151329523807e9 819926.8691406249 … 4.472903950222221e9 -15688.352240668402; -3.8185738483809524e9 -1.7123166484375e6 … -6.2063795222222224e7 119336.55902777778]), NFFTPlan with 284 sampling points for an input array of size(32, 32) and an output array of size(284,) with dims 1:2, ComplexF64[0.9304270810105981 - 0.1755525366839422im, 0.9999999857224657 - 6.37119066797087e-5im, 0.6968566352375541 - 0.5221791050707502im, 0.9304860627318433 + 0.17550202040945986im, 0.8028821861479044 + 0.3950138433443315im, 0.696909133455242 + 0.522149977275377im, 0.9999999859666728 - 6.454986502401626e-5im, 0.690100908057055 - 0.5259755617707915im, 0.9317353097306306 - 0.1646343527093735im, 0.9317937328057632 + 0.1645830409592267im  …  0.8346090489457091 + 0.1795789988845693im, 0.834610800886739 + 0.17958571977595708im, 0.9181402138641586 + 0.09206131734480663im, 0.9120864418553256 + 0.2942607208902332im, 0.9120838428437319 + 0.29426742718619836im, 0.9559925804493471 + 0.09847981391058878im, 0.8881131676530709 + 0.289296040490046im, 0.9924358388500316 + 0.0028007013218525043im, 0.8881114735852544 + 0.28930284794286826im, 0.9999999999233905 + 7.3999334332420275e-6im], BSplinePulse{3}(), [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]), gcache = JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}(sparse([1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [1, 1, 1, 3, 4, 4, 5, 5, 5, 7, 8, 8, 9, 9, 9, 11, 12, 12, 13, 13, 13, 15, 16, 16, 17, 17, 17, 19, 20, 20, 22, 23, 23, 24, 24, 24, 26, 27, 27, 28, 28, 28, 30, 31, 31, 32, 32, 32, 34, 35, 35, 36, 36, 36, 36, 38, 39, 39, 40, 40, 40, 41, 41, 41, 41, 43, 44, 44, 45, 45, 45, 46, 46, 46, 46, 48, 49, 49, 50, 50, 50, 51, 51, 51, 51, 51, 52, 53, 53, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 57, 57, 57, 58, 58, 59, 59, 59, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 64, 64, 64, 64, 64, 64, 65, 65, 66, 66, 66, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 71, 71, 71, 71, 71, 71, 72, 72, 73, 73, 73, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 78, 78, 78, 78, 78, 79, 79, 80, 80, 80, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 86, 86, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 95, 95, 96, 96, 98, 98, 98, 99, 100, 100, 100, 100, 100, 101, 101, 102, 102, 102, 104, 104, 104, 104, 105, 106, 106, 106, 106, 106, 107, 107, 108, 108, 108, 110, 110, 110, 110, 111, 112, 112, 112, 112, 112, 113, 113, 114, 114, 114, 116, 116, 116, 116, 117, 118, 118, 118, 118, 118, 119, 119, 120, 120, 120, 122, 122, 122, 122, 123, 124, 124, 124, 124, 124, 125, 125, 126, 126, 126, 128, 128, 128, 128, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [2, 2, 2, 3, 3, 4, 6, 6, 6, 7, 7, 8, 10, 10, 10, 11, 11, 12, 14, 14, 14, 15, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 22, 25, 25, 25, 26, 26, 27, 29, 29, 29, 30, 30, 31, 33, 33, 33, 34, 34, 35, 37, 37, 37, 37, 38, 38, 38, 39, 39, 40, 42, 42, 42, 42, 43, 43, 43, 44, 44, 45, 47, 47, 47, 47, 48, 48, 48, 49, 49, 50, 52, 52, 52, 52, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 56, 58, 58, 58, 58, 59, 59, 59, 60, 60, 60, 60, 60, 60, 61, 61, 62, 63, 63, 63, 63, 63, 65, 65, 65, 65, 66, 66, 66, 67, 67, 67, 67, 67, 67, 68, 68, 69, 70, 70, 70, 70, 70, 72, 72, 72, 72, 73, 73, 73, 74, 74, 74, 74, 74, 74, 75, 75, 76, 77, 77, 77, 77, 77, 79, 79, 79, 80, 80, 81, 81, 81, 81, 81, 82, 83, 83, 83, 83, 85, 85, 85, 86, 86, 87, 87, 87, 87, 88, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 96, 96, 97, 97, 97, 97, 98, 99, 99, 99, 101, 101, 101, 102, 102, 103, 103, 103, 103, 103, 104, 105, 105, 105, 105, 107, 107, 107, 108, 108, 109, 109, 109, 109, 109, 110, 111, 111, 111, 111, 113, 113, 113, 114, 114, 115, 115, 115, 115, 115, 116, 117, 117, 117, 117, 119, 119, 119, 120, 120, 121, 121, 121, 121, 121, 122, 123, 123, 123, 123, 125, 125, 125, 126, 126, 127, 127, 127, 127, 127, 128, 129, 129, 129, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AA, :AP, :LM, :PV, :AA, :AP, :LM, :PV, :AA, :AP  …  :AZ, :JC, :LM, :SM, :AA, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AA), (0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AA), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AA), (1.516666665673256, :AP)  …  (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AA), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), Fill(NoReference(), 25)), gcachep = JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}(Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([4, 5, 6, 10, 11, 12, 16, 17, 18, 22  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [2, 3, 3, 5, 6, 6, 8, 9, 9, 11, 12, 12, 14, 15, 15, 16, 17, 17, 19, 20, 20, 22, 23, 23, 25, 26, 26, 28, 29, 29, 30, 30, 30, 32, 33, 33, 34, 34, 34, 36, 37, 37, 38, 38, 38, 39, 40, 40, 42, 42, 42, 43, 43, 43, 43, 44, 44, 45, 45, 45, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 50, 50, 51, 51, 51, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 56, 56, 57, 57, 57, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 62, 62, 63, 63, 63, 65, 65, 65, 65, 66, 67, 68, 68, 70, 70, 70, 71, 71, 72, 72, 72, 74, 74, 74, 74, 75, 76, 76, 78, 78, 78, 79, 80, 80, 81, 81, 81, 83, 83, 83, 83, 84, 85, 85, 86, 86, 86, 88, 88, 88, 88, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 96, 96, 96, 98, 98, 98, 98, 99, 100, 100, 101, 101, 101, 103, 103, 103, 103, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [1, 1, 1, 2, 2, 3, 4, 4, 4, 5, 5, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 11, 11, 12, 13, 13, 13, 14, 14, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 21, 22, 22, 23, 24, 24, 24, 25, 25, 26, 27, 27, 27, 27, 28, 28, 28, 29, 29, 30, 31, 31, 31, 31, 32, 32, 32, 33, 33, 34, 35, 35, 35, 35, 36, 36, 36, 37, 37, 38, 39, 39, 39, 39, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 43, 44, 44, 44, 44, 45, 45, 45, 46, 46, 46, 46, 46, 46, 47, 47, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 52, 52, 52, 52, 52, 52, 53, 53, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 58, 58, 58, 58, 58, 58, 59, 59, 60, 61, 61, 61, 61, 61, 62, 62, 62, 63, 63, 64, 64, 64, 64, 64, 65, 66, 66, 66, 66, 67, 67, 67, 68, 68, 69, 69, 69, 69, 70, 71, 71, 71, 72, 72, 73, 73, 73, 73, 73, 74, 75, 75, 75, 75, 76, 76, 77, 77, 77, 77, 78, 79, 79, 79, 80, 80, 80, 81, 81, 82, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 85, 85, 86, 86, 87, 87, 87, 87, 87, 88, 89, 89, 89, 89, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 95, 95, 95, 96, 96, 97, 97, 97, 97, 97, 98, 99, 99, 99, 99, 100, 100, 100, 101, 101, 102, 102, 102, 102, 102, 103, 104, 104, 104, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AP, :LM, :PV, :AP, :LM, :PV, :AP, :LM, :PV, :AP  …  :AP, :AZ, :JC, :LM, :SM, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256, 1.516666665673256, 1.816666603088379  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AP), (1.516666665673256, :LM), (1.516666665673256, :PV), (1.816666603088379, :AP)  …  (7.7166666984558105, :AP), (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), SingleReference{FixedSeg{ComplexF64}}[SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AP, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))  …  SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))]))

We will also fix the total flux to be the observed value 1.1. This is because total flux is degenerate with a global shift in the gain amplitudes making the problem degenerate. To fix this we use the observed total flux as our value.

Moving onto our prior, we first focus on the instrument model priors. Each station requires its own prior on both the amplitudes and phases. For the amplitudes we assume that the gains are apriori well calibrated around unit gains (or 0 log gain amplitudes) which corresponds to no instrument corruption. The gain dispersion is then set to 10% for all stations except LMT, representing that we expect 10% deviations from scan-to-scan. For LMT we let the prior expand to 100% due to the known pointing issues LMT had in 2017.

using Distributions
using DistributionsAD
distamp = station_tuple(dvis, Normal(0.0, 0.1); LM = Normal(1.0))
(AA = Distributions.Normal{Float64}(μ=0.0, σ=0.1), AP = Distributions.Normal{Float64}(μ=0.0, σ=0.1), AZ = Distributions.Normal{Float64}(μ=0.0, σ=0.1), JC = Distributions.Normal{Float64}(μ=0.0, σ=0.1), LM = Distributions.Normal{Float64}(μ=1.0, σ=1.0), PV = Distributions.Normal{Float64}(μ=0.0, σ=0.1), SM = Distributions.Normal{Float64}(μ=0.0, σ=0.1))

For the phases, as mentioned above, we will use a segmented gain prior. This means that rather than the parameters being directly the gains, we fit the first gain for each site, and then the other parameters are the segmented gains compared to the previous time. To model this we break the gain phase prior into two parts. The first is the prior for the first observing timestamp of each site, distphase0, and the second is the prior for segmented gain ϵₜ from time i to i+1, given by distphase. For the EHT, we are dealing with pre-2*rand(rng, ndim) .- 1.5calibrated data, so often, the gain phase jumps from scan to scan are minor. As such, we can put a more informative prior on distphase.

Warning

We use AA (ALMA) as a reference station so we do not have to specify a gain prior for it.

distphase = station_tuple(dvis, DiagonalVonMises(0.0, inv(π^2)))
(AA = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), AP = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), AZ = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), JC = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), LM = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), PV = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), SM = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688))

In addition we want a reasonable guess for what the resolution of our image should be. For radio astronomy this is given by roughly the longest baseline in the image. To put this into pixel space we then divide by the pixel size.

beam = beamsize(dvis)
rat = (beam/(step(grid.X)))
5.279832635689346

To make the Gaussian Markov random field efficient we first precompute a bunch of quantities that allow us to scale things linearly with the number of image pixels. This drastically improves the usual N^3 scaling you get from usual Gaussian Processes.

crcache = MarkovRandomFieldCache(meanpr)
VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}(sparse([1, 2, 32, 33, 993, 1, 2, 3, 34, 994  …  31, 991, 1022, 1023, 1024, 32, 992, 993, 1023, 1024], [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 18, 18, 18, 18, 18, 19, 19, 19, 19, 19, 20, 20, 20, 20, 20, 21, 21, 21, 21, 21, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 24, 24, 24, 24, 24, 25, 25, 25, 25, 25, 26, 26, 26, 26, 26, 27, 27, 27, 27, 27, 28, 28, 28, 28, 28, 29, 29, 29, 29, 29, 30, 30, 30, 30, 30, 31, 31, 31, 31, 31, 32, 32, 32, 32, 32, 33, 33, 33, 33, 33, 34, 34, 34, 34, 34, 35, 35, 35, 35, 35, 36, 36, 36, 36, 36, 37, 37, 37, 37, 37, 38, 38, 38, 38, 38, 39, 39, 39, 39, 39, 40, 40, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 42, 42, 42, 43, 43, 43, 43, 43, 44, 44, 44, 44, 44, 45, 45, 45, 45, 45, 46, 46, 46, 46, 46, 47, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 50, 51, 51, 51, 51, 51, 52, 52, 52, 52, 52, 53, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 56, 57, 57, 57, 57, 57, 58, 58, 58, 58, 58, 59, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 63, 63, 63, 63, 64, 64, 64, 64, 64, 65, 65, 65, 65, 65, 66, 66, 66, 66, 66, 67, 67, 67, 67, 67, 68, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 70, 70, 70, 70, 71, 71, 71, 71, 71, 72, 72, 72, 72, 72, 73, 73, 73, 73, 73, 74, 74, 74, 74, 74, 75, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 77, 77, 77, 77, 78, 78, 78, 78, 78, 79, 79, 79, 79, 79, 80, 80, 80, 80, 80, 81, 81, 81, 81, 81, 82, 82, 82, 82, 82, 83, 83, 83, 83, 83, 84, 84, 84, 84, 84, 85, 85, 85, 85, 85, 86, 86, 86, 86, 86, 87, 87, 87, 87, 87, 88, 88, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 90, 90, 90, 91, 91, 91, 91, 91, 92, 92, 92, 92, 92, 93, 93, 93, 93, 93, 94, 94, 94, 94, 94, 95, 95, 95, 95, 95, 96, 96, 96, 96, 96, 97, 97, 97, 97, 97, 98, 98, 98, 98, 98, 99, 99, 99, 99, 99, 100, 100, 100, 100, 100, 101, 101, 101, 101, 101, 102, 102, 102, 102, 102, 103, 103, 103, 103, 103, 104, 104, 104, 104, 104, 105, 105, 105, 105, 105, 106, 106, 106, 106, 106, 107, 107, 107, 107, 107, 108, 108, 108, 108, 108, 109, 109, 109, 109, 109, 110, 110, 110, 110, 110, 111, 111, 111, 111, 111, 112, 112, 112, 112, 112, 113, 113, 113, 113, 113, 114, 114, 114, 114, 114, 115, 115, 115, 115, 115, 116, 116, 116, 116, 116, 117, 117, 117, 117, 117, 118, 118, 118, 118, 118, 119, 119, 119, 119, 119, 120, 120, 120, 120, 120, 121, 121, 121, 121, 121, 122, 122, 122, 122, 122, 123, 123, 123, 123, 123, 124, 124, 124, 124, 124, 125, 125, 125, 125, 125, 126, 126, 126, 126, 126, 127, 127, 127, 127, 127, 128, 128, 128, 128, 128, 129, 129, 129, 129, 129, 130, 130, 130, 130, 130, 131, 131, 131, 131, 131, 132, 132, 132, 132, 132, 133, 133, 133, 133, 133, 134, 134, 134, 134, 134, 135, 135, 135, 135, 135, 136, 136, 136, 136, 136, 137, 137, 137, 137, 137, 138, 138, 138, 138, 138, 139, 139, 139, 139, 139, 140, 140, 140, 140, 140, 141, 141, 141, 141, 141, 142, 142, 142, 142, 142, 143, 143, 143, 143, 143, 144, 144, 144, 144, 144, 145, 145, 145, 145, 145, 146, 146, 146, 146, 146, 147, 147, 147, 147, 147, 148, 148, 148, 148, 148, 149, 149, 149, 149, 149, 150, 150, 150, 150, 150, 151, 151, 151, 151, 151, 152, 152, 152, 152, 152, 153, 153, 153, 153, 153, 154, 154, 154, 154, 154, 155, 155, 155, 155, 155, 156, 156, 156, 156, 156, 157, 157, 157, 157, 157, 158, 158, 158, 158, 158, 159, 159, 159, 159, 159, 160, 160, 160, 160, 160, 161, 161, 161, 161, 161, 162, 162, 162, 162, 162, 163, 163, 163, 163, 163, 164, 164, 164, 164, 164, 165, 165, 165, 165, 165, 166, 166, 166, 166, 166, 167, 167, 167, 167, 167, 168, 168, 168, 168, 168, 169, 169, 169, 169, 169, 170, 170, 170, 170, 170, 171, 171, 171, 171, 171, 172, 172, 172, 172, 172, 173, 173, 173, 173, 173, 174, 174, 174, 174, 174, 175, 175, 175, 175, 175, 176, 176, 176, 176, 176, 177, 177, 177, 177, 177, 178, 178, 178, 178, 178, 179, 179, 179, 179, 179, 180, 180, 180, 180, 180, 181, 181, 181, 181, 181, 182, 182, 182, 182, 182, 183, 183, 183, 183, 183, 184, 184, 184, 184, 184, 185, 185, 185, 185, 185, 186, 186, 186, 186, 186, 187, 187, 187, 187, 187, 188, 188, 188, 188, 188, 189, 189, 189, 189, 189, 190, 190, 190, 190, 190, 191, 191, 191, 191, 191, 192, 192, 192, 192, 192, 193, 193, 193, 193, 193, 194, 194, 194, 194, 194, 195, 195, 195, 195, 195, 196, 196, 196, 196, 196, 197, 197, 197, 197, 197, 198, 198, 198, 198, 198, 199, 199, 199, 199, 199, 200, 200, 200, 200, 200, 201, 201, 201, 201, 201, 202, 202, 202, 202, 202, 203, 203, 203, 203, 203, 204, 204, 204, 204, 204, 205, 205, 205, 205, 205, 206, 206, 206, 206, 206, 207, 207, 207, 207, 207, 208, 208, 208, 208, 208, 209, 209, 209, 209, 209, 210, 210, 210, 210, 210, 211, 211, 211, 211, 211, 212, 212, 212, 212, 212, 213, 213, 213, 213, 213, 214, 214, 214, 214, 214, 215, 215, 215, 215, 215, 216, 216, 216, 216, 216, 217, 217, 217, 217, 217, 218, 218, 218, 218, 218, 219, 219, 219, 219, 219, 220, 220, 220, 220, 220, 221, 221, 221, 221, 221, 222, 222, 222, 222, 222, 223, 223, 223, 223, 223, 224, 224, 224, 224, 224, 225, 225, 225, 225, 225, 226, 226, 226, 226, 226, 227, 227, 227, 227, 227, 228, 228, 228, 228, 228, 229, 229, 229, 229, 229, 230, 230, 230, 230, 230, 231, 231, 231, 231, 231, 232, 232, 232, 232, 232, 233, 233, 233, 233, 233, 234, 234, 234, 234, 234, 235, 235, 235, 235, 235, 236, 236, 236, 236, 236, 237, 237, 237, 237, 237, 238, 238, 238, 238, 238, 239, 239, 239, 239, 239, 240, 240, 240, 240, 240, 241, 241, 241, 241, 241, 242, 242, 242, 242, 242, 243, 243, 243, 243, 243, 244, 244, 244, 244, 244, 245, 245, 245, 245, 245, 246, 246, 246, 246, 246, 247, 247, 247, 247, 247, 248, 248, 248, 248, 248, 249, 249, 249, 249, 249, 250, 250, 250, 250, 250, 251, 251, 251, 251, 251, 252, 252, 252, 252, 252, 253, 253, 253, 253, 253, 254, 254, 254, 254, 254, 255, 255, 255, 255, 255, 256, 256, 256, 256, 256, 257, 257, 257, 257, 257, 258, 258, 258, 258, 258, 259, 259, 259, 259, 259, 260, 260, 260, 260, 260, 261, 261, 261, 261, 261, 262, 262, 262, 262, 262, 263, 263, 263, 263, 263, 264, 264, 264, 264, 264, 265, 265, 265, 265, 265, 266, 266, 266, 266, 266, 267, 267, 267, 267, 267, 268, 268, 268, 268, 268, 269, 269, 269, 269, 269, 270, 270, 270, 270, 270, 271, 271, 271, 271, 271, 272, 272, 272, 272, 272, 273, 273, 273, 273, 273, 274, 274, 274, 274, 274, 275, 275, 275, 275, 275, 276, 276, 276, 276, 276, 277, 277, 277, 277, 277, 278, 278, 278, 278, 278, 279, 279, 279, 279, 279, 280, 280, 280, 280, 280, 281, 281, 281, 281, 281, 282, 282, 282, 282, 282, 283, 283, 283, 283, 283, 284, 284, 284, 284, 284, 285, 285, 285, 285, 285, 286, 286, 286, 286, 286, 287, 287, 287, 287, 287, 288, 288, 288, 288, 288, 289, 289, 289, 289, 289, 290, 290, 290, 290, 290, 291, 291, 291, 291, 291, 292, 292, 292, 292, 292, 293, 293, 293, 293, 293, 294, 294, 294, 294, 294, 295, 295, 295, 295, 295, 296, 296, 296, 296, 296, 297, 297, 297, 297, 297, 298, 298, 298, 298, 298, 299, 299, 299, 299, 299, 300, 300, 300, 300, 300, 301, 301, 301, 301, 301, 302, 302, 302, 302, 302, 303, 303, 303, 303, 303, 304, 304, 304, 304, 304, 305, 305, 305, 305, 305, 306, 306, 306, 306, 306, 307, 307, 307, 307, 307, 308, 308, 308, 308, 308, 309, 309, 309, 309, 309, 310, 310, 310, 310, 310, 311, 311, 311, 311, 311, 312, 312, 312, 312, 312, 313, 313, 313, 313, 313, 314, 314, 314, 314, 314, 315, 315, 315, 315, 315, 316, 316, 316, 316, 316, 317, 317, 317, 317, 317, 318, 318, 318, 318, 318, 319, 319, 319, 319, 319, 320, 320, 320, 320, 320, 321, 321, 321, 321, 321, 322, 322, 322, 322, 322, 323, 323, 323, 323, 323, 324, 324, 324, 324, 324, 325, 325, 325, 325, 325, 326, 326, 326, 326, 326, 327, 327, 327, 327, 327, 328, 328, 328, 328, 328, 329, 329, 329, 329, 329, 330, 330, 330, 330, 330, 331, 331, 331, 331, 331, 332, 332, 332, 332, 332, 333, 333, 333, 333, 333, 334, 334, 334, 334, 334, 335, 335, 335, 335, 335, 336, 336, 336, 336, 336, 337, 337, 337, 337, 337, 338, 338, 338, 338, 338, 339, 339, 339, 339, 339, 340, 340, 340, 340, 340, 341, 341, 341, 341, 341, 342, 342, 342, 342, 342, 343, 343, 343, 343, 343, 344, 344, 344, 344, 344, 345, 345, 345, 345, 345, 346, 346, 346, 346, 346, 347, 347, 347, 347, 347, 348, 348, 348, 348, 348, 349, 349, 349, 349, 349, 350, 350, 350, 350, 350, 351, 351, 351, 351, 351, 352, 352, 352, 352, 352, 353, 353, 353, 353, 353, 354, 354, 354, 354, 354, 355, 355, 355, 355, 355, 356, 356, 356, 356, 356, 357, 357, 357, 357, 357, 358, 358, 358, 358, 358, 359, 359, 359, 359, 359, 360, 360, 360, 360, 360, 361, 361, 361, 361, 361, 362, 362, 362, 362, 362, 363, 363, 363, 363, 363, 364, 364, 364, 364, 364, 365, 365, 365, 365, 365, 366, 366, 366, 366, 366, 367, 367, 367, 367, 367, 368, 368, 368, 368, 368, 369, 369, 369, 369, 369, 370, 370, 370, 370, 370, 371, 371, 371, 371, 371, 372, 372, 372, 372, 372, 373, 373, 373, 373, 373, 374, 374, 374, 374, 374, 375, 375, 375, 375, 375, 376, 376, 376, 376, 376, 377, 377, 377, 377, 377, 378, 378, 378, 378, 378, 379, 379, 379, 379, 379, 380, 380, 380, 380, 380, 381, 381, 381, 381, 381, 382, 382, 382, 382, 382, 383, 383, 383, 383, 383, 384, 384, 384, 384, 384, 385, 385, 385, 385, 385, 386, 386, 386, 386, 386, 387, 387, 387, 387, 387, 388, 388, 388, 388, 388, 389, 389, 389, 389, 389, 390, 390, 390, 390, 390, 391, 391, 391, 391, 391, 392, 392, 392, 392, 392, 393, 393, 393, 393, 393, 394, 394, 394, 394, 394, 395, 395, 395, 395, 395, 396, 396, 396, 396, 396, 397, 397, 397, 397, 397, 398, 398, 398, 398, 398, 399, 399, 399, 399, 399, 400, 400, 400, 400, 400, 401, 401, 401, 401, 401, 402, 402, 402, 402, 402, 403, 403, 403, 403, 403, 404, 404, 404, 404, 404, 405, 405, 405, 405, 405, 406, 406, 406, 406, 406, 407, 407, 407, 407, 407, 408, 408, 408, 408, 408, 409, 409, 409, 409, 409, 410, 410, 410, 410, 410, 411, 411, 411, 411, 411, 412, 412, 412, 412, 412, 413, 413, 413, 413, 413, 414, 414, 414, 414, 414, 415, 415, 415, 415, 415, 416, 416, 416, 416, 416, 417, 417, 417, 417, 417, 418, 418, 418, 418, 418, 419, 419, 419, 419, 419, 420, 420, 420, 420, 420, 421, 421, 421, 421, 421, 422, 422, 422, 422, 422, 423, 423, 423, 423, 423, 424, 424, 424, 424, 424, 425, 425, 425, 425, 425, 426, 426, 426, 426, 426, 427, 427, 427, 427, 427, 428, 428, 428, 428, 428, 429, 429, 429, 429, 429, 430, 430, 430, 430, 430, 431, 431, 431, 431, 431, 432, 432, 432, 432, 432, 433, 433, 433, 433, 433, 434, 434, 434, 434, 434, 435, 435, 435, 435, 435, 436, 436, 436, 436, 436, 437, 437, 437, 437, 437, 438, 438, 438, 438, 438, 439, 439, 439, 439, 439, 440, 440, 440, 440, 440, 441, 441, 441, 441, 441, 442, 442, 442, 442, 442, 443, 443, 443, 443, 443, 444, 444, 444, 444, 444, 445, 445, 445, 445, 445, 446, 446, 446, 446, 446, 447, 447, 447, 447, 447, 448, 448, 448, 448, 448, 449, 449, 449, 449, 449, 450, 450, 450, 450, 450, 451, 451, 451, 451, 451, 452, 452, 452, 452, 452, 453, 453, 453, 453, 453, 454, 454, 454, 454, 454, 455, 455, 455, 455, 455, 456, 456, 456, 456, 456, 457, 457, 457, 457, 457, 458, 458, 458, 458, 458, 459, 459, 459, 459, 459, 460, 460, 460, 460, 460, 461, 461, 461, 461, 461, 462, 462, 462, 462, 462, 463, 463, 463, 463, 463, 464, 464, 464, 464, 464, 465, 465, 465, 465, 465, 466, 466, 466, 466, 466, 467, 467, 467, 467, 467, 468, 468, 468, 468, 468, 469, 469, 469, 469, 469, 470, 470, 470, 470, 470, 471, 471, 471, 471, 471, 472, 472, 472, 472, 472, 473, 473, 473, 473, 473, 474, 474, 474, 474, 474, 475, 475, 475, 475, 475, 476, 476, 476, 476, 476, 477, 477, 477, 477, 477, 478, 478, 478, 478, 478, 479, 479, 479, 479, 479, 480, 480, 480, 480, 480, 481, 481, 481, 481, 481, 482, 482, 482, 482, 482, 483, 483, 483, 483, 483, 484, 484, 484, 484, 484, 485, 485, 485, 485, 485, 486, 486, 486, 486, 486, 487, 487, 487, 487, 487, 488, 488, 488, 488, 488, 489, 489, 489, 489, 489, 490, 490, 490, 490, 490, 491, 491, 491, 491, 491, 492, 492, 492, 492, 492, 493, 493, 493, 493, 493, 494, 494, 494, 494, 494, 495, 495, 495, 495, 495, 496, 496, 496, 496, 496, 497, 497, 497, 497, 497, 498, 498, 498, 498, 498, 499, 499, 499, 499, 499, 500, 500, 500, 500, 500, 501, 501, 501, 501, 501, 502, 502, 502, 502, 502, 503, 503, 503, 503, 503, 504, 504, 504, 504, 504, 505, 505, 505, 505, 505, 506, 506, 506, 506, 506, 507, 507, 507, 507, 507, 508, 508, 508, 508, 508, 509, 509, 509, 509, 509, 510, 510, 510, 510, 510, 511, 511, 511, 511, 511, 512, 512, 512, 512, 512, 513, 513, 513, 513, 513, 514, 514, 514, 514, 514, 515, 515, 515, 515, 515, 516, 516, 516, 516, 516, 517, 517, 517, 517, 517, 518, 518, 518, 518, 518, 519, 519, 519, 519, 519, 520, 520, 520, 520, 520, 521, 521, 521, 521, 521, 522, 522, 522, 522, 522, 523, 523, 523, 523, 523, 524, 524, 524, 524, 524, 525, 525, 525, 525, 525, 526, 526, 526, 526, 526, 527, 527, 527, 527, 527, 528, 528, 528, 528, 528, 529, 529, 529, 529, 529, 530, 530, 530, 530, 530, 531, 531, 531, 531, 531, 532, 532, 532, 532, 532, 533, 533, 533, 533, 533, 534, 534, 534, 534, 534, 535, 535, 535, 535, 535, 536, 536, 536, 536, 536, 537, 537, 537, 537, 537, 538, 538, 538, 538, 538, 539, 539, 539, 539, 539, 540, 540, 540, 540, 540, 541, 541, 541, 541, 541, 542, 542, 542, 542, 542, 543, 543, 543, 543, 543, 544, 544, 544, 544, 544, 545, 545, 545, 545, 545, 546, 546, 546, 546, 546, 547, 547, 547, 547, 547, 548, 548, 548, 548, 548, 549, 549, 549, 549, 549, 550, 550, 550, 550, 550, 551, 551, 551, 551, 551, 552, 552, 552, 552, 552, 553, 553, 553, 553, 553, 554, 554, 554, 554, 554, 555, 555, 555, 555, 555, 556, 556, 556, 556, 556, 557, 557, 557, 557, 557, 558, 558, 558, 558, 558, 559, 559, 559, 559, 559, 560, 560, 560, 560, 560, 561, 561, 561, 561, 561, 562, 562, 562, 562, 562, 563, 563, 563, 563, 563, 564, 564, 564, 564, 564, 565, 565, 565, 565, 565, 566, 566, 566, 566, 566, 567, 567, 567, 567, 567, 568, 568, 568, 568, 568, 569, 569, 569, 569, 569, 570, 570, 570, 570, 570, 571, 571, 571, 571, 571, 572, 572, 572, 572, 572, 573, 573, 573, 573, 573, 574, 574, 574, 574, 574, 575, 575, 575, 575, 575, 576, 576, 576, 576, 576, 577, 577, 577, 577, 577, 578, 578, 578, 578, 578, 579, 579, 579, 579, 579, 580, 580, 580, 580, 580, 581, 581, 581, 581, 581, 582, 582, 582, 582, 582, 583, 583, 583, 583, 583, 584, 584, 584, 584, 584, 585, 585, 585, 585, 585, 586, 586, 586, 586, 586, 587, 587, 587, 587, 587, 588, 588, 588, 588, 588, 589, 589, 589, 589, 589, 590, 590, 590, 590, 590, 591, 591, 591, 591, 591, 592, 592, 592, 592, 592, 593, 593, 593, 593, 593, 594, 594, 594, 594, 594, 595, 595, 595, 595, 595, 596, 596, 596, 596, 596, 597, 597, 597, 597, 597, 598, 598, 598, 598, 598, 599, 599, 599, 599, 599, 600, 600, 600, 600, 600, 601, 601, 601, 601, 601, 602, 602, 602, 602, 602, 603, 603, 603, 603, 603, 604, 604, 604, 604, 604, 605, 605, 605, 605, 605, 606, 606, 606, 606, 606, 607, 607, 607, 607, 607, 608, 608, 608, 608, 608, 609, 609, 609, 609, 609, 610, 610, 610, 610, 610, 611, 611, 611, 611, 611, 612, 612, 612, 612, 612, 613, 613, 613, 613, 613, 614, 614, 614, 614, 614, 615, 615, 615, 615, 615, 616, 616, 616, 616, 616, 617, 617, 617, 617, 617, 618, 618, 618, 618, 618, 619, 619, 619, 619, 619, 620, 620, 620, 620, 620, 621, 621, 621, 621, 621, 622, 622, 622, 622, 622, 623, 623, 623, 623, 623, 624, 624, 624, 624, 624, 625, 625, 625, 625, 625, 626, 626, 626, 626, 626, 627, 627, 627, 627, 627, 628, 628, 628, 628, 628, 629, 629, 629, 629, 629, 630, 630, 630, 630, 630, 631, 631, 631, 631, 631, 632, 632, 632, 632, 632, 633, 633, 633, 633, 633, 634, 634, 634, 634, 634, 635, 635, 635, 635, 635, 636, 636, 636, 636, 636, 637, 637, 637, 637, 637, 638, 638, 638, 638, 638, 639, 639, 639, 639, 639, 640, 640, 640, 640, 640, 641, 641, 641, 641, 641, 642, 642, 642, 642, 642, 643, 643, 643, 643, 643, 644, 644, 644, 644, 644, 645, 645, 645, 645, 645, 646, 646, 646, 646, 646, 647, 647, 647, 647, 647, 648, 648, 648, 648, 648, 649, 649, 649, 649, 649, 650, 650, 650, 650, 650, 651, 651, 651, 651, 651, 652, 652, 652, 652, 652, 653, 653, 653, 653, 653, 654, 654, 654, 654, 654, 655, 655, 655, 655, 655, 656, 656, 656, 656, 656, 657, 657, 657, 657, 657, 658, 658, 658, 658, 658, 659, 659, 659, 659, 659, 660, 660, 660, 660, 660, 661, 661, 661, 661, 661, 662, 662, 662, 662, 662, 663, 663, 663, 663, 663, 664, 664, 664, 664, 664, 665, 665, 665, 665, 665, 666, 666, 666, 666, 666, 667, 667, 667, 667, 667, 668, 668, 668, 668, 668, 669, 669, 669, 669, 669, 670, 670, 670, 670, 670, 671, 671, 671, 671, 671, 672, 672, 672, 672, 672, 673, 673, 673, 673, 673, 674, 674, 674, 674, 674, 675, 675, 675, 675, 675, 676, 676, 676, 676, 676, 677, 677, 677, 677, 677, 678, 678, 678, 678, 678, 679, 679, 679, 679, 679, 680, 680, 680, 680, 680, 681, 681, 681, 681, 681, 682, 682, 682, 682, 682, 683, 683, 683, 683, 683, 684, 684, 684, 684, 684, 685, 685, 685, 685, 685, 686, 686, 686, 686, 686, 687, 687, 687, 687, 687, 688, 688, 688, 688, 688, 689, 689, 689, 689, 689, 690, 690, 690, 690, 690, 691, 691, 691, 691, 691, 692, 692, 692, 692, 692, 693, 693, 693, 693, 693, 694, 694, 694, 694, 694, 695, 695, 695, 695, 695, 696, 696, 696, 696, 696, 697, 697, 697, 697, 697, 698, 698, 698, 698, 698, 699, 699, 699, 699, 699, 700, 700, 700, 700, 700, 701, 701, 701, 701, 701, 702, 702, 702, 702, 702, 703, 703, 703, 703, 703, 704, 704, 704, 704, 704, 705, 705, 705, 705, 705, 706, 706, 706, 706, 706, 707, 707, 707, 707, 707, 708, 708, 708, 708, 708, 709, 709, 709, 709, 709, 710, 710, 710, 710, 710, 711, 711, 711, 711, 711, 712, 712, 712, 712, 712, 713, 713, 713, 713, 713, 714, 714, 714, 714, 714, 715, 715, 715, 715, 715, 716, 716, 716, 716, 716, 717, 717, 717, 717, 717, 718, 718, 718, 718, 718, 719, 719, 719, 719, 719, 720, 720, 720, 720, 720, 721, 721, 721, 721, 721, 722, 722, 722, 722, 722, 723, 723, 723, 723, 723, 724, 724, 724, 724, 724, 725, 725, 725, 725, 725, 726, 726, 726, 726, 726, 727, 727, 727, 727, 727, 728, 728, 728, 728, 728, 729, 729, 729, 729, 729, 730, 730, 730, 730, 730, 731, 731, 731, 731, 731, 732, 732, 732, 732, 732, 733, 733, 733, 733, 733, 734, 734, 734, 734, 734, 735, 735, 735, 735, 735, 736, 736, 736, 736, 736, 737, 737, 737, 737, 737, 738, 738, 738, 738, 738, 739, 739, 739, 739, 739, 740, 740, 740, 740, 740, 741, 741, 741, 741, 741, 742, 742, 742, 742, 742, 743, 743, 743, 743, 743, 744, 744, 744, 744, 744, 745, 745, 745, 745, 745, 746, 746, 746, 746, 746, 747, 747, 747, 747, 747, 748, 748, 748, 748, 748, 749, 749, 749, 749, 749, 750, 750, 750, 750, 750, 751, 751, 751, 751, 751, 752, 752, 752, 752, 752, 753, 753, 753, 753, 753, 754, 754, 754, 754, 754, 755, 755, 755, 755, 755, 756, 756, 756, 756, 756, 757, 757, 757, 757, 757, 758, 758, 758, 758, 758, 759, 759, 759, 759, 759, 760, 760, 760, 760, 760, 761, 761, 761, 761, 761, 762, 762, 762, 762, 762, 763, 763, 763, 763, 763, 764, 764, 764, 764, 764, 765, 765, 765, 765, 765, 766, 766, 766, 766, 766, 767, 767, 767, 767, 767, 768, 768, 768, 768, 768, 769, 769, 769, 769, 769, 770, 770, 770, 770, 770, 771, 771, 771, 771, 771, 772, 772, 772, 772, 772, 773, 773, 773, 773, 773, 774, 774, 774, 774, 774, 775, 775, 775, 775, 775, 776, 776, 776, 776, 776, 777, 777, 777, 777, 777, 778, 778, 778, 778, 778, 779, 779, 779, 779, 779, 780, 780, 780, 780, 780, 781, 781, 781, 781, 781, 782, 782, 782, 782, 782, 783, 783, 783, 783, 783, 784, 784, 784, 784, 784, 785, 785, 785, 785, 785, 786, 786, 786, 786, 786, 787, 787, 787, 787, 787, 788, 788, 788, 788, 788, 789, 789, 789, 789, 789, 790, 790, 790, 790, 790, 791, 791, 791, 791, 791, 792, 792, 792, 792, 792, 793, 793, 793, 793, 793, 794, 794, 794, 794, 794, 795, 795, 795, 795, 795, 796, 796, 796, 796, 796, 797, 797, 797, 797, 797, 798, 798, 798, 798, 798, 799, 799, 799, 799, 799, 800, 800, 800, 800, 800, 801, 801, 801, 801, 801, 802, 802, 802, 802, 802, 803, 803, 803, 803, 803, 804, 804, 804, 804, 804, 805, 805, 805, 805, 805, 806, 806, 806, 806, 806, 807, 807, 807, 807, 807, 808, 808, 808, 808, 808, 809, 809, 809, 809, 809, 810, 810, 810, 810, 810, 811, 811, 811, 811, 811, 812, 812, 812, 812, 812, 813, 813, 813, 813, 813, 814, 814, 814, 814, 814, 815, 815, 815, 815, 815, 816, 816, 816, 816, 816, 817, 817, 817, 817, 817, 818, 818, 818, 818, 818, 819, 819, 819, 819, 819, 820, 820, 820, 820, 820, 821, 821, 821, 821, 821, 822, 822, 822, 822, 822, 823, 823, 823, 823, 823, 824, 824, 824, 824, 824, 825, 825, 825, 825, 825, 826, 826, 826, 826, 826, 827, 827, 827, 827, 827, 828, 828, 828, 828, 828, 829, 829, 829, 829, 829, 830, 830, 830, 830, 830, 831, 831, 831, 831, 831, 832, 832, 832, 832, 832, 833, 833, 833, 833, 833, 834, 834, 834, 834, 834, 835, 835, 835, 835, 835, 836, 836, 836, 836, 836, 837, 837, 837, 837, 837, 838, 838, 838, 838, 838, 839, 839, 839, 839, 839, 840, 840, 840, 840, 840, 841, 841, 841, 841, 841, 842, 842, 842, 842, 842, 843, 843, 843, 843, 843, 844, 844, 844, 844, 844, 845, 845, 845, 845, 845, 846, 846, 846, 846, 846, 847, 847, 847, 847, 847, 848, 848, 848, 848, 848, 849, 849, 849, 849, 849, 850, 850, 850, 850, 850, 851, 851, 851, 851, 851, 852, 852, 852, 852, 852, 853, 853, 853, 853, 853, 854, 854, 854, 854, 854, 855, 855, 855, 855, 855, 856, 856, 856, 856, 856, 857, 857, 857, 857, 857, 858, 858, 858, 858, 858, 859, 859, 859, 859, 859, 860, 860, 860, 860, 860, 861, 861, 861, 861, 861, 862, 862, 862, 862, 862, 863, 863, 863, 863, 863, 864, 864, 864, 864, 864, 865, 865, 865, 865, 865, 866, 866, 866, 866, 866, 867, 867, 867, 867, 867, 868, 868, 868, 868, 868, 869, 869, 869, 869, 869, 870, 870, 870, 870, 870, 871, 871, 871, 871, 871, 872, 872, 872, 872, 872, 873, 873, 873, 873, 873, 874, 874, 874, 874, 874, 875, 875, 875, 875, 875, 876, 876, 876, 876, 876, 877, 877, 877, 877, 877, 878, 878, 878, 878, 878, 879, 879, 879, 879, 879, 880, 880, 880, 880, 880, 881, 881, 881, 881, 881, 882, 882, 882, 882, 882, 883, 883, 883, 883, 883, 884, 884, 884, 884, 884, 885, 885, 885, 885, 885, 886, 886, 886, 886, 886, 887, 887, 887, 887, 887, 888, 888, 888, 888, 888, 889, 889, 889, 889, 889, 890, 890, 890, 890, 890, 891, 891, 891, 891, 891, 892, 892, 892, 892, 892, 893, 893, 893, 893, 893, 894, 894, 894, 894, 894, 895, 895, 895, 895, 895, 896, 896, 896, 896, 896, 897, 897, 897, 897, 897, 898, 898, 898, 898, 898, 899, 899, 899, 899, 899, 900, 900, 900, 900, 900, 901, 901, 901, 901, 901, 902, 902, 902, 902, 902, 903, 903, 903, 903, 903, 904, 904, 904, 904, 904, 905, 905, 905, 905, 905, 906, 906, 906, 906, 906, 907, 907, 907, 907, 907, 908, 908, 908, 908, 908, 909, 909, 909, 909, 909, 910, 910, 910, 910, 910, 911, 911, 911, 911, 911, 912, 912, 912, 912, 912, 913, 913, 913, 913, 913, 914, 914, 914, 914, 914, 915, 915, 915, 915, 915, 916, 916, 916, 916, 916, 917, 917, 917, 917, 917, 918, 918, 918, 918, 918, 919, 919, 919, 919, 919, 920, 920, 920, 920, 920, 921, 921, 921, 921, 921, 922, 922, 922, 922, 922, 923, 923, 923, 923, 923, 924, 924, 924, 924, 924, 925, 925, 925, 925, 925, 926, 926, 926, 926, 926, 927, 927, 927, 927, 927, 928, 928, 928, 928, 928, 929, 929, 929, 929, 929, 930, 930, 930, 930, 930, 931, 931, 931, 931, 931, 932, 932, 932, 932, 932, 933, 933, 933, 933, 933, 934, 934, 934, 934, 934, 935, 935, 935, 935, 935, 936, 936, 936, 936, 936, 937, 937, 937, 937, 937, 938, 938, 938, 938, 938, 939, 939, 939, 939, 939, 940, 940, 940, 940, 940, 941, 941, 941, 941, 941, 942, 942, 942, 942, 942, 943, 943, 943, 943, 943, 944, 944, 944, 944, 944, 945, 945, 945, 945, 945, 946, 946, 946, 946, 946, 947, 947, 947, 947, 947, 948, 948, 948, 948, 948, 949, 949, 949, 949, 949, 950, 950, 950, 950, 950, 951, 951, 951, 951, 951, 952, 952, 952, 952, 952, 953, 953, 953, 953, 953, 954, 954, 954, 954, 954, 955, 955, 955, 955, 955, 956, 956, 956, 956, 956, 957, 957, 957, 957, 957, 958, 958, 958, 958, 958, 959, 959, 959, 959, 959, 960, 960, 960, 960, 960, 961, 961, 961, 961, 961, 962, 962, 962, 962, 962, 963, 963, 963, 963, 963, 964, 964, 964, 964, 964, 965, 965, 965, 965, 965, 966, 966, 966, 966, 966, 967, 967, 967, 967, 967, 968, 968, 968, 968, 968, 969, 969, 969, 969, 969, 970, 970, 970, 970, 970, 971, 971, 971, 971, 971, 972, 972, 972, 972, 972, 973, 973, 973, 973, 973, 974, 974, 974, 974, 974, 975, 975, 975, 975, 975, 976, 976, 976, 976, 976, 977, 977, 977, 977, 977, 978, 978, 978, 978, 978, 979, 979, 979, 979, 979, 980, 980, 980, 980, 980, 981, 981, 981, 981, 981, 982, 982, 982, 982, 982, 983, 983, 983, 983, 983, 984, 984, 984, 984, 984, 985, 985, 985, 985, 985, 986, 986, 986, 986, 986, 987, 987, 987, 987, 987, 988, 988, 988, 988, 988, 989, 989, 989, 989, 989, 990, 990, 990, 990, 990, 991, 991, 991, 991, 991, 992, 992, 992, 992, 992, 993, 993, 993, 993, 993, 994, 994, 994, 994, 994, 995, 995, 995, 995, 995, 996, 996, 996, 996, 996, 997, 997, 997, 997, 997, 998, 998, 998, 998, 998, 999, 999, 999, 999, 999, 1000, 1000, 1000, 1000, 1000, 1001, 1001, 1001, 1001, 1001, 1002, 1002, 1002, 1002, 1002, 1003, 1003, 1003, 1003, 1003, 1004, 1004, 1004, 1004, 1004, 1005, 1005, 1005, 1005, 1005, 1006, 1006, 1006, 1006, 1006, 1007, 1007, 1007, 1007, 1007, 1008, 1008, 1008, 1008, 1008, 1009, 1009, 1009, 1009, 1009, 1010, 1010, 1010, 1010, 1010, 1011, 1011, 1011, 1011, 1011, 1012, 1012, 1012, 1012, 1012, 1013, 1013, 1013, 1013, 1013, 1014, 1014, 1014, 1014, 1014, 1015, 1015, 1015, 1015, 1015, 1016, 1016, 1016, 1016, 1016, 1017, 1017, 1017, 1017, 1017, 1018, 1018, 1018, 1018, 1018, 1019, 1019, 1019, 1019, 1019, 1020, 1020, 1020, 1020, 1020, 1021, 1021, 1021, 1021, 1021, 1022, 1022, 1022, 1022, 1022, 1023, 1023, 1023, 1023, 1023, 1024, 1024, 1024, 1024, 1024], [4.0, -1.0, -1.0, -1.0, -1.0, -1.0, 4.0, -1.0, -1.0, -1.0  …  -1.0, -1.0, -1.0, 4.0, -1.0, -1.0, -1.0, -1.0, -1.0, 4.0], 1024, 1024), [1.0 0.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 1.0], [0.0 0.03842943919353914 … 0.15224093497742697 0.03842943919353936; 0.03842943919353914 0.07685887838707828 … 0.1906703741709661 0.0768588783870785; … ; 0.15224093497742697 0.1906703741709661 … 0.30448186995485393 0.19067037417096633; 0.03842943919353914 0.07685887838707828 … 0.1906703741709661 0.0768588783870785])

One of the benefits of the Bayesian approach is that we can fit for the hyperparameters of our prior/regularizers unlike traditional RML appraoches. To construct this heirarchical prior we will first make a map that takes in our regularizer hyperparameters and returns the image prior given those hyperparameters.

fmap = let crcache=crcache
    x->GaussMarkovRandomField(x, 1.0, crcache)
end
#1 (generic function with 1 method)

Now we can finally form our image prior. For this we use a heirarchical prior where the inverse correlation length is given by a Half-Normal distribution whose peak is at zero and standard deviation is 0.1/rat where recall rat is the beam size per pixel. For the variance of the random field we use another half normal prior with standard deviation 0.1. The reason we use the half-normal priors is to prefer "simple" structures. Gaussian Markov random fields are extremly flexible models, and to prevent overfitting it is common to use priors that penalize complexity. Therefore, we want to use priors that enforce similarity to our mean image. If the data wants more complexity then it will drive us away from the prior.

cprior = HierarchicalPrior(fmap, InverseGamma(1.0, -log(0.01*rat)))
VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}(
priormap: #1
hyperprior: Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.3399885194615155)
θ: 2.9412757865584154
)

)

We can now form our model parameter priors. Like our other imaging examples, we use a Dirichlet prior for our image pixels. For the log gain amplitudes, we use the CalPrior which automatically constructs the prior for the given jones cache gcache.

prior = NamedDist(
         c = cprior,
         fg = Uniform(0.0, 1.0),
         σimg = truncated(Normal(0.0, 0.1); lower = 0.0),
         lgamp = CalPrior(distamp, gcache),
         gphase = CalPrior(distphase, gcachep),
        )
VLBIImagePriors.NamedDist{(:c, :fg, :σimg, :lgamp, :gphase), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Uniform{Float64}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}(
dists: (VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}(
priormap: #1
hyperprior: Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.3399885194615155)
θ: 2.9412757865584154
)

)
, Distributions.Uniform{Float64}(a=0.0, b=1.0), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.1); lower=0.0), CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}(
dists: DiagNormal(
dim: 129
μ: [0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0]
Σ: [0.010000000000000002 0.0 … 0.0 0.0; 0.0 0.010000000000000002 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 0.010000000000000002]
)

jcache: JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}(sparse([1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [1, 1, 1, 3, 4, 4, 5, 5, 5, 7, 8, 8, 9, 9, 9, 11, 12, 12, 13, 13, 13, 15, 16, 16, 17, 17, 17, 19, 20, 20, 22, 23, 23, 24, 24, 24, 26, 27, 27, 28, 28, 28, 30, 31, 31, 32, 32, 32, 34, 35, 35, 36, 36, 36, 36, 38, 39, 39, 40, 40, 40, 41, 41, 41, 41, 43, 44, 44, 45, 45, 45, 46, 46, 46, 46, 48, 49, 49, 50, 50, 50, 51, 51, 51, 51, 51, 52, 53, 53, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 57, 57, 57, 58, 58, 59, 59, 59, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 64, 64, 64, 64, 64, 64, 65, 65, 66, 66, 66, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 71, 71, 71, 71, 71, 71, 72, 72, 73, 73, 73, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 78, 78, 78, 78, 78, 79, 79, 80, 80, 80, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 86, 86, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 95, 95, 96, 96, 98, 98, 98, 99, 100, 100, 100, 100, 100, 101, 101, 102, 102, 102, 104, 104, 104, 104, 105, 106, 106, 106, 106, 106, 107, 107, 108, 108, 108, 110, 110, 110, 110, 111, 112, 112, 112, 112, 112, 113, 113, 114, 114, 114, 116, 116, 116, 116, 117, 118, 118, 118, 118, 118, 119, 119, 120, 120, 120, 122, 122, 122, 122, 123, 124, 124, 124, 124, 124, 125, 125, 126, 126, 126, 128, 128, 128, 128, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [2, 2, 2, 3, 3, 4, 6, 6, 6, 7, 7, 8, 10, 10, 10, 11, 11, 12, 14, 14, 14, 15, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 22, 25, 25, 25, 26, 26, 27, 29, 29, 29, 30, 30, 31, 33, 33, 33, 34, 34, 35, 37, 37, 37, 37, 38, 38, 38, 39, 39, 40, 42, 42, 42, 42, 43, 43, 43, 44, 44, 45, 47, 47, 47, 47, 48, 48, 48, 49, 49, 50, 52, 52, 52, 52, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 56, 58, 58, 58, 58, 59, 59, 59, 60, 60, 60, 60, 60, 60, 61, 61, 62, 63, 63, 63, 63, 63, 65, 65, 65, 65, 66, 66, 66, 67, 67, 67, 67, 67, 67, 68, 68, 69, 70, 70, 70, 70, 70, 72, 72, 72, 72, 73, 73, 73, 74, 74, 74, 74, 74, 74, 75, 75, 76, 77, 77, 77, 77, 77, 79, 79, 79, 80, 80, 81, 81, 81, 81, 81, 82, 83, 83, 83, 83, 85, 85, 85, 86, 86, 87, 87, 87, 87, 88, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 96, 96, 97, 97, 97, 97, 98, 99, 99, 99, 101, 101, 101, 102, 102, 103, 103, 103, 103, 103, 104, 105, 105, 105, 105, 107, 107, 107, 108, 108, 109, 109, 109, 109, 109, 110, 111, 111, 111, 111, 113, 113, 113, 114, 114, 115, 115, 115, 115, 115, 116, 117, 117, 117, 117, 119, 119, 119, 120, 120, 121, 121, 121, 121, 121, 122, 123, 123, 123, 123, 125, 125, 125, 126, 126, 127, 127, 127, 127, 127, 128, 129, 129, 129, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AA, :AP, :LM, :PV, :AA, :AP, :LM, :PV, :AA, :AP  …  :AZ, :JC, :LM, :SM, :AA, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AA), (0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AA), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AA), (1.516666665673256, :AP)  …  (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AA), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), Fill(NoReference(), 25))
)
, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}(
dists: VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}(
μ: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
κ: [0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778  …  0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778]
lnorm: -180.86855628209554
)

jcache: JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}(Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([4, 5, 6, 10, 11, 12, 16, 17, 18, 22  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [2, 3, 3, 5, 6, 6, 8, 9, 9, 11, 12, 12, 14, 15, 15, 16, 17, 17, 19, 20, 20, 22, 23, 23, 25, 26, 26, 28, 29, 29, 30, 30, 30, 32, 33, 33, 34, 34, 34, 36, 37, 37, 38, 38, 38, 39, 40, 40, 42, 42, 42, 43, 43, 43, 43, 44, 44, 45, 45, 45, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 50, 50, 51, 51, 51, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 56, 56, 57, 57, 57, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 62, 62, 63, 63, 63, 65, 65, 65, 65, 66, 67, 68, 68, 70, 70, 70, 71, 71, 72, 72, 72, 74, 74, 74, 74, 75, 76, 76, 78, 78, 78, 79, 80, 80, 81, 81, 81, 83, 83, 83, 83, 84, 85, 85, 86, 86, 86, 88, 88, 88, 88, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 96, 96, 96, 98, 98, 98, 98, 99, 100, 100, 101, 101, 101, 103, 103, 103, 103, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [1, 1, 1, 2, 2, 3, 4, 4, 4, 5, 5, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 11, 11, 12, 13, 13, 13, 14, 14, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 21, 22, 22, 23, 24, 24, 24, 25, 25, 26, 27, 27, 27, 27, 28, 28, 28, 29, 29, 30, 31, 31, 31, 31, 32, 32, 32, 33, 33, 34, 35, 35, 35, 35, 36, 36, 36, 37, 37, 38, 39, 39, 39, 39, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 43, 44, 44, 44, 44, 45, 45, 45, 46, 46, 46, 46, 46, 46, 47, 47, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 52, 52, 52, 52, 52, 52, 53, 53, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 58, 58, 58, 58, 58, 58, 59, 59, 60, 61, 61, 61, 61, 61, 62, 62, 62, 63, 63, 64, 64, 64, 64, 64, 65, 66, 66, 66, 66, 67, 67, 67, 68, 68, 69, 69, 69, 69, 70, 71, 71, 71, 72, 72, 73, 73, 73, 73, 73, 74, 75, 75, 75, 75, 76, 76, 77, 77, 77, 77, 78, 79, 79, 79, 80, 80, 80, 81, 81, 82, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 85, 85, 86, 86, 87, 87, 87, 87, 87, 88, 89, 89, 89, 89, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 95, 95, 95, 96, 96, 97, 97, 97, 97, 97, 98, 99, 99, 99, 99, 100, 100, 100, 101, 101, 102, 102, 102, 102, 102, 103, 104, 104, 104, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AP, :LM, :PV, :AP, :LM, :PV, :AP, :LM, :PV, :AP  …  :AP, :AZ, :JC, :LM, :SM, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256, 1.516666665673256, 1.816666603088379  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AP), (1.516666665673256, :LM), (1.516666665673256, :PV), (1.816666603088379, :AP)  …  (7.7166666984558105, :AP), (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), SingleReference{FixedSeg{ComplexF64}}[SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AP, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))  …  SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))])
)
)
)

Putting it all together we form our likelihood and posterior objects for optimization and sampling.

lklhd = RadioLikelihood(sky, instrument, dvis; skymeta=metadata, instrumentmeta=metadata)
post = Posterior(lklhd, prior)
Posterior{RadioLikelihood{Comrade.ModelMetadata{typeof(Main.sky), NamedTuple{(:ftot, :K, :meanpr, :grid, :cache, :gcache, :gcachep), Tuple{Float64, VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}, Matrix{Float64}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}, VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}, Comrade.ModelMetadata{typeof(Main.instrument), NamedTuple{(:ftot, :K, :meanpr, :grid, :cache, :gcache, :gcachep), Tuple{Float64, VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}, Matrix{Float64}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}, VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}, Tuple{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}}, Tuple{Comrade.ConditionedLikelihood{Comrade.var"#26#27"{Float64, Vector{Float64}}, Vector{ComplexF64}}}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, NamedTuple{(:U, :V, :T, :F), NTuple{4, Vector{Float64}}}}, VLBIImagePriors.NamedDist{(:c, :fg, :σimg, :lgamp, :gphase), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Uniform{Float64}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}}(RadioLikelihood
	Number of data products: 1
, VLBIImagePriors.NamedDist{(:c, :fg, :σimg, :lgamp, :gphase), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Uniform{Float64}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}(
dists: (VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}(
priormap: #1
hyperprior: Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.3399885194615155)
θ: 2.9412757865584154
)

)
, Distributions.Uniform{Float64}(a=0.0, b=1.0), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.1); lower=0.0), CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}(
dists: DiagNormal(
dim: 129
μ: [0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0]
Σ: [0.010000000000000002 0.0 … 0.0 0.0; 0.0 0.010000000000000002 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 0.010000000000000002]
)

jcache: JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}(sparse([1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [1, 1, 1, 3, 4, 4, 5, 5, 5, 7, 8, 8, 9, 9, 9, 11, 12, 12, 13, 13, 13, 15, 16, 16, 17, 17, 17, 19, 20, 20, 22, 23, 23, 24, 24, 24, 26, 27, 27, 28, 28, 28, 30, 31, 31, 32, 32, 32, 34, 35, 35, 36, 36, 36, 36, 38, 39, 39, 40, 40, 40, 41, 41, 41, 41, 43, 44, 44, 45, 45, 45, 46, 46, 46, 46, 48, 49, 49, 50, 50, 50, 51, 51, 51, 51, 51, 52, 53, 53, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 57, 57, 57, 58, 58, 59, 59, 59, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 64, 64, 64, 64, 64, 64, 65, 65, 66, 66, 66, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 71, 71, 71, 71, 71, 71, 72, 72, 73, 73, 73, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 78, 78, 78, 78, 78, 79, 79, 80, 80, 80, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 86, 86, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 95, 95, 96, 96, 98, 98, 98, 99, 100, 100, 100, 100, 100, 101, 101, 102, 102, 102, 104, 104, 104, 104, 105, 106, 106, 106, 106, 106, 107, 107, 108, 108, 108, 110, 110, 110, 110, 111, 112, 112, 112, 112, 112, 113, 113, 114, 114, 114, 116, 116, 116, 116, 117, 118, 118, 118, 118, 118, 119, 119, 120, 120, 120, 122, 122, 122, 122, 123, 124, 124, 124, 124, 124, 125, 125, 126, 126, 126, 128, 128, 128, 128, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [2, 2, 2, 3, 3, 4, 6, 6, 6, 7, 7, 8, 10, 10, 10, 11, 11, 12, 14, 14, 14, 15, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 22, 25, 25, 25, 26, 26, 27, 29, 29, 29, 30, 30, 31, 33, 33, 33, 34, 34, 35, 37, 37, 37, 37, 38, 38, 38, 39, 39, 40, 42, 42, 42, 42, 43, 43, 43, 44, 44, 45, 47, 47, 47, 47, 48, 48, 48, 49, 49, 50, 52, 52, 52, 52, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 56, 58, 58, 58, 58, 59, 59, 59, 60, 60, 60, 60, 60, 60, 61, 61, 62, 63, 63, 63, 63, 63, 65, 65, 65, 65, 66, 66, 66, 67, 67, 67, 67, 67, 67, 68, 68, 69, 70, 70, 70, 70, 70, 72, 72, 72, 72, 73, 73, 73, 74, 74, 74, 74, 74, 74, 75, 75, 76, 77, 77, 77, 77, 77, 79, 79, 79, 80, 80, 81, 81, 81, 81, 81, 82, 83, 83, 83, 83, 85, 85, 85, 86, 86, 87, 87, 87, 87, 88, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 96, 96, 97, 97, 97, 97, 98, 99, 99, 99, 101, 101, 101, 102, 102, 103, 103, 103, 103, 103, 104, 105, 105, 105, 105, 107, 107, 107, 108, 108, 109, 109, 109, 109, 109, 110, 111, 111, 111, 111, 113, 113, 113, 114, 114, 115, 115, 115, 115, 115, 116, 117, 117, 117, 117, 119, 119, 119, 120, 120, 121, 121, 121, 121, 121, 122, 123, 123, 123, 123, 125, 125, 125, 126, 126, 127, 127, 127, 127, 127, 128, 129, 129, 129, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AA, :AP, :LM, :PV, :AA, :AP, :LM, :PV, :AA, :AP  …  :AZ, :JC, :LM, :SM, :AA, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AA), (0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AA), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AA), (1.516666665673256, :AP)  …  (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AA), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), Fill(NoReference(), 25))
)
, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}(
dists: VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}(
μ: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
κ: [0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778  …  0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778]
lnorm: -180.86855628209554
)

jcache: JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}(Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([4, 5, 6, 10, 11, 12, 16, 17, 18, 22  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [2, 3, 3, 5, 6, 6, 8, 9, 9, 11, 12, 12, 14, 15, 15, 16, 17, 17, 19, 20, 20, 22, 23, 23, 25, 26, 26, 28, 29, 29, 30, 30, 30, 32, 33, 33, 34, 34, 34, 36, 37, 37, 38, 38, 38, 39, 40, 40, 42, 42, 42, 43, 43, 43, 43, 44, 44, 45, 45, 45, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 50, 50, 51, 51, 51, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 56, 56, 57, 57, 57, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 62, 62, 63, 63, 63, 65, 65, 65, 65, 66, 67, 68, 68, 70, 70, 70, 71, 71, 72, 72, 72, 74, 74, 74, 74, 75, 76, 76, 78, 78, 78, 79, 80, 80, 81, 81, 81, 83, 83, 83, 83, 84, 85, 85, 86, 86, 86, 88, 88, 88, 88, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 96, 96, 96, 98, 98, 98, 98, 99, 100, 100, 101, 101, 101, 103, 103, 103, 103, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [1, 1, 1, 2, 2, 3, 4, 4, 4, 5, 5, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 11, 11, 12, 13, 13, 13, 14, 14, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 21, 22, 22, 23, 24, 24, 24, 25, 25, 26, 27, 27, 27, 27, 28, 28, 28, 29, 29, 30, 31, 31, 31, 31, 32, 32, 32, 33, 33, 34, 35, 35, 35, 35, 36, 36, 36, 37, 37, 38, 39, 39, 39, 39, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 43, 44, 44, 44, 44, 45, 45, 45, 46, 46, 46, 46, 46, 46, 47, 47, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 52, 52, 52, 52, 52, 52, 53, 53, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 58, 58, 58, 58, 58, 58, 59, 59, 60, 61, 61, 61, 61, 61, 62, 62, 62, 63, 63, 64, 64, 64, 64, 64, 65, 66, 66, 66, 66, 67, 67, 67, 68, 68, 69, 69, 69, 69, 70, 71, 71, 71, 72, 72, 73, 73, 73, 73, 73, 74, 75, 75, 75, 75, 76, 76, 77, 77, 77, 77, 78, 79, 79, 79, 80, 80, 80, 81, 81, 82, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 85, 85, 86, 86, 87, 87, 87, 87, 87, 88, 89, 89, 89, 89, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 95, 95, 95, 96, 96, 97, 97, 97, 97, 97, 98, 99, 99, 99, 99, 100, 100, 100, 101, 101, 102, 102, 102, 102, 102, 103, 104, 104, 104, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AP, :LM, :PV, :AP, :LM, :PV, :AP, :LM, :PV, :AP  …  :AP, :AZ, :JC, :LM, :SM, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256, 1.516666665673256, 1.816666603088379  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AP), (1.516666665673256, :LM), (1.516666665673256, :PV), (1.816666603088379, :AP)  …  (7.7166666984558105, :AP), (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), SingleReference{FixedSeg{ComplexF64}}[SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AP, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))  …  SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))])
)
)
)
)

Reconstructing the Image and Instrument Effects

To sample from this posterior, it is convenient to move from our constrained parameter space to an unconstrained one (i.e., the support of the transformed posterior is (-∞, ∞)). This is done using the asflat function.

tpost = asflat(post)
ndim = dimension(tpost)
1364

Our Posterior and TransformedPosterior objects satisfy the LogDensityProblems interface. This allows us to easily switch between different AD backends and many of Julia's statistical inference packages use this interface as well.

using LogDensityProblemsAD
using Zygote
gtpost = ADgradient(Val(:Zygote), tpost)
x0 = randn(rng, ndim)
LogDensityProblemsAD.logdensity_and_gradient(gtpost, x0)
(-1.0889677227706838e8, [3.9714863653464385, -41.82940853200137, 69.59557261663501, 4.569805366219421, 10.789295400079078, 240.91028613995184, 39.14697994202466, 42.2142944805375, 50.171331063617025, 69.99422943458495  …  3666.7681059308334, 11816.75911472515, 1587.1490272240208, -343.3074914007036, 6760.359870612455, 349.8400750594293, 17534.646730431534, -7352.639076547881, -181.3718982383253, 410.6393309543062])

We can now also find the dimension of our posterior or the number of parameters we are going to sample.

Warning

This can often be different from what you would expect. This is especially true when using angular variables where we often artificially increase the dimension of the parameter space to make sampling easier.

To initialize our sampler we will use optimize using LBFGS

using ComradeOptimization
using OptimizationOptimJL
f = OptimizationFunction(tpost, Optimization.AutoZygote())
prob = Optimization.OptimizationProblem(f, prior_sample(rng, tpost), nothing)
ℓ = logdensityof(tpost)
sol = solve(prob, LBFGS(), maxiters=4_000, g_tol=1e-1);

Now transform back to parameter space

xopt = transform(tpost, sol.u)
(c = (params = [-6.8795954061674e-6 -7.857487634196908e-6 … -1.3035633065842715e-6 -2.129336117648613e-6; -1.1099349857742423e-5 -1.8778000230961743e-5 … -1.050227717163778e-5 -5.493128246328554e-6; … ; 8.607840915016035e-6 4.865994508196023e-6 … 7.083591345017068e-6 2.2022154708638837e-6; 2.1423414747582233e-6 2.0005035322394562e-7 … 7.068683632421048e-6 1.4368280706136605e-6], hyperparams = 0.05624068712337777), fg = 0.005847386813875843, σimg = 3.106827059593349, lgamp = [-0.008994995924661119, 0.008974868355736848, -0.8872560828253759, 0.13657360396825097, 0.014431082559460275, -0.01395683722815999, -0.3160012148065026, 0.24217916493353425, 0.005625711920769715, -0.005245328088601519  …  -0.2913372838812374, -0.0006601217815102756, -1.0958333041747068, -0.003033343152992906, 0.026366114788167937, -0.025657185444572624, -0.2803354713864622, 0.0005445606468547087, -1.154881241510689, -0.004427900326664654], gphase = [-0.8142694298905171, 3.135755867015113, 2.348870625945773, -0.8711505416498114, -3.0495456438519066, 2.1548697081081594, -0.928549610360143, -2.9630003959697078, 1.9653167711040551, -0.9773398807580612  …  -0.5846227667564549, -0.5034836443985393, -0.5312567371376251, -0.698818717558769, 1.2401582641783033, -0.5238623533119972, -0.6160567263019271, -0.5001654780284682, -0.9106415275813774, 1.288858981207776])
Warning

Fitting gains tends to be very difficult, meaning that optimization can take a lot longer. The upside is that we usually get nicer images.

First we will evaluate our fit by plotting the residuals

using Plots
residual(vlbimodel(post, xopt), dvis)
Example block output

These look reasonable, although there may be some minor overfitting. This could be improved in a few ways, but that is beyond the goal of this quick tutorial. Plotting the image, we see that we have a much cleaner version of the closure-only image from Imaging a Black Hole using only Closure Quantities.

import WGLMakie as CM
img = intensitymap(skymodel(post, xopt), fovx, fovy, 128, 128)
CM.image(img, axis=(xreversed=true, aspect=1, title="MAP Image"), colormap=:afmhot)

Because we also fit the instrument model, we can inspect their parameters. To do this, Comrade provides a caltable function that converts the flattened gain parameters to a tabular format based on the time and its segmentation.

gt = Comrade.caltable(gcachep, xopt.gphase)
plot(gt, layout=(3,3), size=(600,500))
Example block output

The gain phases are pretty random, although much of this is due to us picking a random reference station for each scan.

Moving onto the gain amplitudes, we see that most of the gain variation is within 10% as expected except LMT, which has massive variations.

gt = Comrade.caltable(gcache, exp.(xopt.lgamp))
plot(gt, layout=(3,3), size=(600,500))
Example block output

To sample from the posterior, we will use HMC, specifically the NUTS algorithm. For information about NUTS, see Michael Betancourt's notes.

Note

For our metric, we use a diagonal matrix due to easier tuning

However, due to the need to sample a large number of gain parameters, constructing the posterior is rather time-consuming. Therefore, for this tutorial, we will only do a quick preliminary run, and any posterior inferences should be appropriately skeptical.

using ComradeAHMC
metric = DiagEuclideanMetric(ndim)
chain, stats = sample(rng, post, AHMC(;metric, autodiff=Val(:Zygote)), 700; nadapts=500, init_params=xopt)
(NamedTuple{(:c, :fg, :σimg, :lgamp, :gphase), Tuple{NamedTuple{(:params, :hyperparams), Tuple{Matrix{Float64}, Float64}}, Float64, Float64, Vector{Float64}, Vector{Float64}}}[(c = (params = [-0.0015773160299475209 0.007446292868229237 … 0.01170766171051302 -0.007005680391699723; 0.007162992124767366 0.013366960848094479 … 0.00015862047659218019 -0.003304582632558772; … ; -0.005644354818011102 0.0044400607394182775 … -0.004532897602210734 0.0016241492621931282; -0.009683625514645177 0.010529549529139115 … -0.0035474376423090346 -0.0033065157410006574], hyperparams = 0.05631323874081243), fg = 0.005912899727762349, σimg = 3.121275496378236, lgamp = [-0.009704333221575092, 0.009128961387364406, -0.8911202231241673, 0.1343924273013784, 0.021070460529475633, -0.020958175206573476, -0.3242547898802818, 0.22950450900385058, 0.004428637732820791, -0.003946341697229905  …  -0.2926454022266335, 0.0019077225329123466, -1.1029708564257206, 0.0015870506943462936, 0.023956200188540184, -0.02355129944396012, -0.2823117712302015, -0.006928285429001772, -1.158773871006694, 0.0030844149785100613], gphase = [-0.8173263812513099, -3.131603168288999, 2.35720951172687, -0.8685770196334339, -3.0614266186628316, 2.1636969398322985, -0.9280506965404098, -2.964578408476249, 1.9686986618955808, -0.978748356367405  …  -0.5869689846395485, -0.5100773657323292, -0.5175244870834038, -0.6991341656956106, 1.2511799139838944, -0.5219043981098191, -0.6137093520635004, -0.502127740876625, -0.9132378881286314, 1.282731485668361]), (c = (params = [-0.0015773160299475209 0.007446292868229237 … 0.01170766171051302 -0.007005680391699723; 0.007162992124767366 0.013366960848094479 … 0.00015862047659218019 -0.003304582632558772; … ; -0.005644354818011102 0.0044400607394182775 … -0.004532897602210734 0.0016241492621931282; -0.009683625514645177 0.010529549529139115 … -0.0035474376423090346 -0.0033065157410006574], hyperparams = 0.05631323874081243), fg = 0.005912899727762349, σimg = 3.121275496378236, lgamp = [-0.009704333221575092, 0.009128961387364406, -0.8911202231241673, 0.1343924273013784, 0.021070460529475633, -0.020958175206573476, -0.3242547898802818, 0.22950450900385058, 0.004428637732820791, -0.003946341697229905  …  -0.2926454022266335, 0.0019077225329123466, -1.1029708564257206, 0.0015870506943462936, 0.023956200188540184, -0.02355129944396012, -0.2823117712302015, -0.006928285429001772, -1.158773871006694, 0.0030844149785100613], gphase = [-0.8173263812513099, -3.131603168288999, 2.35720951172687, -0.8685770196334339, -3.0614266186628316, 2.1636969398322985, -0.9280506965404098, -2.964578408476249, 1.9686986618955808, -0.978748356367405  …  -0.5869689846395485, -0.5100773657323292, -0.5175244870834038, -0.6991341656956106, 1.2511799139838944, -0.5219043981098191, -0.6137093520635004, -0.502127740876625, -0.9132378881286314, 1.282731485668361]), (c = (params = [-0.11457967227027577 -0.013822420133175589 … -0.17025753382156028 -0.12933332738996875; 0.012523959557493353 -0.024003476535508963 … -0.05000212515001226 -0.02398232175518964; … ; -0.08892953959226617 0.0010620177757260584 … -0.10518362236095523 -0.09351517710582011; -0.10918112801995306 -0.035652865253860445 … -0.05833917484366804 0.0618934144015261], hyperparams = 0.11931318583874885), fg = 0.004725678525854248, σimg = 2.5467442189392386, lgamp = [0.002024561635571309, -0.005029803304642936, -0.8851483595587311, 0.11074991268958634, -0.0014572296746859243, 0.003396600357660719, -0.31146240068203324, 0.24817858636949708, 0.00021673320862473434, -0.001012495491291064  …  -0.2990951286552984, -0.00658979276030965, -1.1150149343572253, -0.01125180657506084, 0.025613327997798675, -0.024442754875115314, -0.2753252978018315, -0.0008167994553909488, -1.1656762931622742, -0.0032732395444089297], gphase = [-0.8187179286982674, 3.1179824689288904, 2.377205796290804, -0.8665205769793897, -3.0457652032129405, 2.182505188560306, -0.9292935666220893, -2.9593118514095478, 1.993240825958968, -0.9775973807893158  …  -0.5811467239901583, -0.5107157894408257, -0.5521958444676569, -0.7003010239707789, 1.219113434160715, -0.5211166738669102, -0.6208381197003475, -0.5321011703822826, -0.904564951828671, 1.2481282382751333]), (c = (params = [0.042226290044198436 0.09725149797895653 … -0.30631255241763616 0.0013302388201500005; -0.025438701614754096 -0.1909063060864143 … -0.16695013660470154 0.04996156316703059; … ; 0.03188545485667665 -0.04066721100295436 … -0.06591494754772642 -0.10996204234572918; -0.1460647243484858 -0.09095851890002524 … 0.010602605662457252 0.11703102679859394], hyperparams = 0.16994500665506407), fg = 0.005801351771989194, σimg = 2.0226995909639855, lgamp = [0.0017182365567845854, -0.0018559657999831418, -0.9892679981015122, 0.0720346954718507, -0.004923938919374527, 0.001079965191310559, -0.3850928196923113, 0.2061771572582628, -0.00877396094206329, 0.0018252512870239003  …  -0.2902895113627469, -0.017580509301092995, -1.1135490951060685, 5.200537268675151e-5, 0.008886677824797705, -0.012788471605683718, -0.2681606357415109, 0.008435510139463477, -1.1682687750546747, -0.010622892048835953], gphase = [-0.8152893147017696, -3.096536222837162, 2.3682562227290607, -0.8670659539213434, -3.0134112586106694, 2.179208092706323, -0.9333562728356278, -2.9524298496283214, 1.9916320707692652, -0.9770063658608801  …  -0.5870481978157686, -0.4921416470582145, -0.5943570922470337, -0.6928778057246447, 1.177891622472667, -0.5274125311971337, -0.6185262171720917, -0.5354580658169265, -0.9049276825093804, 1.247549849298824]), (c = (params = [-0.11209495123522492 -0.13680625876839453 … 0.3671142929124708 0.04423613697359076; 0.03860791420559574 0.23127191710440323 … 0.16350584779847188 -0.11125518301061443; … ; 0.021187304812146472 -0.0013678763699381902 … 0.02454780200283731 0.16304719362499637; 0.2285187295969955 0.08049196999061543 … -0.07216392451959786 -0.0892323476044677], hyperparams = 0.24102034272439346), fg = 0.02034050162086382, σimg = 1.5832757105949082, lgamp = [0.015719094601165003, -0.01753464609084777, -1.114775554157464, -0.04715066777597266, -0.012847415488167042, 0.0076316502364407995, -0.49675174016064544, 0.11131782265827443, -0.023363656514966503, 0.024740869737771874  …  -0.2568159124243884, 0.014565336714222665, -1.1309346943541856, -0.015223080849162476, -0.005728170831801161, 0.010127558957407332, -0.25411942334417387, -0.010741130154318653, -1.2026091368239993, -0.0007751161411654956], gphase = [-0.8171152074104616, -3.1240460115584985, 2.2339928665318736, -0.8746817100624815, -3.0379973577083, 2.0258003862995597, -0.9292607142602416, -2.958502029072243, 1.8378189617460883, -0.9807307124676002  …  -0.5868634977889817, -0.5124104084382244, -0.4140856614385768, -0.6904916832336491, 1.3555943470506666, -0.5313655928351193, -0.617098494407559, -0.37167247546411136, -0.8910344786574561, 1.4206601016583813]), (c = (params = [0.11166361767732208 -0.03244775132790386 … -0.4358284211316636 -0.012366089430065806; -0.2878645502762428 -0.38001280527234665 … -0.28842041431515447 -0.14785942741321648; … ; -0.050761050864243956 0.047105970921332556 … 0.02547698726211179 -0.17150392006712203; -0.3551409122458712 -0.15636785280563148 … -0.10850334062881782 -0.051172115787584675], hyperparams = 0.26333986989406344), fg = 0.036874277013884814, σimg = 1.5453687971443606, lgamp = [-0.028279113852210328, 0.025845326990737612, -1.0611738499847945, -0.02770164385633856, 0.008471670296246094, -0.0059597859405027855, -0.45964166572489645, 0.07250963719210772, -0.01589808466385068, 0.013043490722830545  …  -0.24942817871767148, -0.016390599065013944, -1.129762772861323, 0.00015360484355660304, 0.001737632079942872, -0.000832851085807573, -0.2623236184421148, 0.016489492466901078, -1.1973659863663242, -0.023842902200206486], gphase = [-0.8111847611276605, 3.1202352833072253, 2.3318363801078896, -0.8656102076811905, -3.0590396037263865, 2.126636405075343, -0.9310212588771293, -2.9862098818861074, 1.9474837333713768, -0.978436008217914  …  -0.5862604242058439, -0.5480740330872846, -0.43551681514747065, -0.7288258219527585, 1.3189644349492813, -0.5287144696596441, -0.636748528318168, -0.42395143597338114, -0.9326816997811078, 1.3762788526620329]), (c = (params = [0.11330324715282303 -0.060090203873232946 … -0.43016186064444023 0.021127409270485308; -0.30491090778822993 -0.4117473783453391 … -0.30307624297828445 -0.144900248178917; … ; -0.0370883066340404 0.009361546652933124 … 0.015028397933505857 -0.17247046276561603; -0.3438825568615835 -0.1765712851541505 … -0.10940002629739955 -0.005877062704222146], hyperparams = 0.26565089455561147), fg = 0.03656923663045306, σimg = 1.5283109027066515, lgamp = [-0.0393021181362009, 0.03551526482769094, -1.007858453901983, -0.006038030788676423, -0.008362324252082096, 0.011760112019503558, -0.4419287398499371, 0.09734699298259396, -0.012289465656894321, 0.013647191811190694  …  -0.26588000513863347, -0.0009451563031676385, -1.1172614638479668, 0.005546922071458119, 0.002176016248855465, -0.0016238650255787448, -0.2626127644956636, 0.01748361613897028, -1.1961016461778908, -0.01808772735866742], gphase = [-0.8096445818084903, 3.137648190087009, 2.3323472452775302, -0.8725193828898784, -3.055050753193334, 2.140990162480258, -0.9319525717923256, -2.9930142796947217, 1.9474259864544274, -0.9776222908933424  …  -0.58795789782554, -0.5516463490699479, -0.4725823509678281, -0.7329239411833897, 1.3220769228338325, -0.5180605350831263, -0.6516475265378823, -0.3974372936454414, -0.9320260422605079, 1.3794058333544466]), (c = (params = [0.2763534370573807 -0.18089356585406163 … -0.19295005901437845 -0.018335535706438994; -0.27101212189551144 -0.3769025794558457 … -0.25835877248795674 -0.16055707823823634; … ; -0.27508336663189964 -0.027641747216536918 … -0.050521173684548976 -0.1910435956761797; -0.40851649290134007 -0.20793885889365207 … -0.23963810095340926 -0.11158321596957728], hyperparams = 0.30284037175449224), fg = 0.03572405534059923, σimg = 1.540376275074309, lgamp = [0.042966323824126164, -0.03778056661393969, -1.1060756980438122, -0.030807939678973408, 0.003959693454522108, -0.009047095921354628, -0.438677727318679, 0.14679849311140802, -0.009304991991825027, 0.006308568545170375  …  -0.26636842523241366, 0.007435300674883907, -1.0943529549906037, 0.002358557268523694, 0.011261881685490045, -0.012036414349390234, -0.26100907673058826, -0.015065531812101445, -1.184040183198659, 0.005204499916103467], gphase = [-0.821361571813222, -3.121370645903419, 2.280792189002498, -0.8716208114633418, -3.026426140458346, 2.099914861386887, -0.9250042235088197, -2.9670494168152164, 1.8956406176517668, -0.9769939688760215  …  -0.5842264577865572, -0.5307467861967992, -0.39863461620186375, -0.717300037465643, 1.3713185459455604, -0.5260543674743842, -0.664684231836143, -0.39246897036811684, -0.9216662714325676, 1.4110346981010202]), (c = (params = [0.18345498709976407 -0.0606077934647376 … 0.103489376594723 0.1600065470859707; -0.31321368502089936 -0.34417967822904855 … -0.41133553845310533 -0.07644746935704598; … ; -0.3862809068719504 0.03972100178736266 … 0.14919949001369137 -0.13082114023661584; -0.3367732549536559 -0.12032297425022409 … -0.2194462751511841 0.12231547629234278], hyperparams = 0.2735315211086454), fg = 0.03207968371617366, σimg = 1.5959141045974603, lgamp = [-0.08106142041411935, 0.07906795914006268, -1.005812351690794, 0.07294927979565385, -0.010017288603332447, 0.010886381000207675, -0.4709403435677224, 0.13872933698279197, -0.010253821408314047, 0.01172228102714113  …  -0.2849085433443794, 0.004098571257175896, -1.1075284128746512, -0.003406258615960989, 0.01484817393387975, -0.015285879616584007, -0.24865609740814434, 0.01857977705808561, -1.1976263979449788, -0.020181890696986417], gphase = [-0.8124623589581144, -3.1121022981832387, 2.3038482471415893, -0.8702624159512045, -3.0383783553090535, 2.107289805005533, -0.9247918680621481, -2.988367197818669, 1.9346525113843633, -0.9778130257038552  …  -0.5862215844461068, -0.5363653173760348, -0.4722672835721622, -0.7272742671101189, 1.314835344519491, -0.5199514150521418, -0.6382219406408387, -0.4016855091490881, -0.9301585339744544, 1.3834985922403344]), (c = (params = [0.14008755567777034 -0.08640894419396165 … -0.16570178618297654 -0.05645708891250884; -0.36218590126150507 -0.35328987252701993 … 0.3109845253310053 0.21089477186800315; … ; 0.19276147543438113 -0.16397594805180302 … -0.11135755922883159 -0.38274576873023075; -0.1473463860393902 0.2655364968679685 … -0.20067119074575557 -0.14268153465917607], hyperparams = 0.42291233098359576), fg = 0.05551364886585114, σimg = 1.5074369925339244, lgamp = [0.00516139152596486, -0.002876771563678985, -1.1385377868484297, -0.04002320327216473, -0.0269289128154537, 0.027524394464817612, -0.5237536462580924, 0.12213351748651974, -0.02531291342274148, 0.02464069725731534  …  -0.23518755715089323, 0.011312030433882356, -1.1382643188544475, -0.012425277777840523, -0.030961870195344472, 0.03270327236640663, -0.2179033423234322, -0.00041199000783345534, -1.1971778672088533, -0.00554771429302845], gphase = [-0.8165084080345296, -3.006285214803882, 2.343405784430073, -0.8693019688032644, -2.929750880755671, 2.1613254242529676, -0.9287553536222989, -2.8812663115717174, 1.9846018307762001, -0.9786318404107097  …  -0.5824142017400127, -0.46741803503086227, -0.34600676465376284, -0.663378499142617, 1.4149354613518592, -0.5252057746098926, -0.5908637446836497, -0.32527366737241015, -0.8677240807467568, 1.4665040640133542])  …  (c = (params = [-0.49181236851645965 -1.011685403723889 … -0.6238671823037727 -0.047898257144772534; 0.13385865250049395 1.073578534390941 … -0.6132584426126267 -0.6775980930552307; … ; 0.15392951551434508 -0.6565643172574549 … -0.5249118808171435 0.8435901448466043; -0.6607421906018751 -0.49496793978132986 … 0.4001136159031692 0.9407868727314227], hyperparams = 44.495083024456136), fg = 0.28063537105346364, σimg = 0.8652575939816205, lgamp = [0.016224966856890476, -0.021470596090791758, -0.8809318364224663, -0.06794318591982496, -0.033286500217615304, 0.02445570690272776, -0.2520428847567521, 0.1256795487576371, -0.029545469056276086, 0.032772131303124036  …  -0.08938547167199734, -0.0033568763466695355, -0.9135885091512423, 0.0005528882766075417, -0.01233348349230561, 0.006582950191574444, -0.1040700237895878, -0.001397331879574564, -0.9387638033123215, -0.00050066048864683], gphase = [-0.811141475524782, -2.2952643171034617, 1.106951219283878, -0.8699349797853374, -2.242624431317048, 0.915468217258222, -0.9278696259549382, -2.1481684280205298, 0.7420075076418021, -0.9783000183211137  …  -0.5838530618983229, 0.47379617524244994, 1.7439102937709714, 0.08604380614256217, -2.7640379129571073, -0.5283344442126272, 0.30789401391253085, 1.7361977542093807, -0.15518270676962148, -2.7534627584846296]), (c = (params = [-0.9272365860440605 -1.126661165223037 … -0.3436378540288673 -0.05099910362384746; 0.3650039628547554 0.3837966498699662 … -0.3098955763188251 -0.7374882381820144; … ; 0.7059254880460316 -0.19404000769555946 … -0.3902951571755544 1.1789449063418251; -0.5664229999357817 -0.28951932165413713 … 0.18046109223691442 0.5497961210695541], hyperparams = 40.18523412251972), fg = 0.2860435147201816, σimg = 0.8997465993831403, lgamp = [0.019122783463051542, -0.017270959625966417, -0.8147285630336826, -0.07242575930801536, -0.023614517347618998, 0.02254149541076703, -0.2205547808243384, 0.11892561225765855, -0.02497110522025951, 0.02514432002992846  …  -0.09292092283055373, -0.012735831322363871, -0.9139246959251646, -0.0038522214253139896, -0.011218493695877064, 0.017148631793104253, -0.10556954704646995, -0.004926409988026401, -0.9666022442665109, 0.0014584955770888686], gphase = [-0.810051337880092, -2.325287502979104, 1.0912599820690772, -0.8710871882952397, -2.2457329821693204, 0.8930302362418585, -0.9281327478025293, -2.165266714630528, 0.7121242123440878, -0.9732537631977157  …  -0.5862210375603776, 0.45421977942261743, 1.7488114973164923, 0.0624054236509602, -2.762706459697255, -0.5282722226555652, 0.3126318443185473, 1.7419631790724321, -0.18209426646142957, -2.7554846752727244]), (c = (params = [0.4811566656956373 0.3624749713630575 … -0.5797184524524794 0.2078903377357365; 0.5738619497072355 0.3565130729010635 … 0.9121015128685035 0.3504824624102584; … ; 1.1643769819488967 0.47628161862549595 … -1.6323329019042279 -0.09762446394682421; -0.3800912176535576 1.6673825351252032 … -0.7775763072830127 -0.7393597005139521], hyperparams = 25.894614544255322), fg = 0.26102665574490147, σimg = 0.9754115440874708, lgamp = [-0.042070608638390475, 0.037913772307171335, -0.8133692814034381, -0.053477093761341146, 0.006377412548316855, -0.016342281855667523, -0.33663761874835435, 0.043521080636755104, 0.0013425083704176788, 0.00013149283455338843  …  -0.09849245823146917, -0.004479158012405493, -0.98249673302851, -0.015127307011833757, 0.007902411343622085, -0.010384230203891526, -0.13393102968157383, 0.0012821853751536635, -1.0326492451866727, -0.007461812229053084], gphase = [-0.8176064064311312, -2.5729878347470003, 0.7797011640053312, -0.8710626557799953, -2.475829576403006, 0.5836772893806323, -0.9287334831854713, -2.3858873870965347, 0.39960546503592176, -0.9750920066010587  …  -0.5830331858270835, 0.2150643186007219, 1.6528688526313782, -0.14187126072978432, -2.867439105083193, -0.5255125842642813, 0.060192106508226845, 1.6278116697590352, -0.38088162086915817, -2.865909748578458]), (c = (params = [0.5844913511260921 0.2951783567885804 … 0.6700396080880833 1.56749913839649; 0.34252984310142404 0.723967738572967 … 0.946155732717254 0.5246897959716865; … ; 1.1152235696135249 -0.37903360049358953 … 1.4232353362076309 1.5145581195728615; -0.5523382773219144 0.04624237729427234 … 1.8610279613971974 1.3270361214229474], hyperparams = 104.52200122113119), fg = 0.290742025665108, σimg = 0.9100408099345091, lgamp = [0.03162549415319128, -0.03284844163853763, -0.9254649082893074, -0.03360439039004699, -0.016331271966292926, 0.01945318908650759, -0.3224775310322246, 0.14401006268370714, -0.02387522383874977, 0.025230753834022145  …  -0.03959846226551538, -0.012636551252886683, -0.9733963794218936, -0.009007118342957006, -0.024861565922390787, 0.025491556035920755, -0.06006264029869964, 0.0007051668697634268, -1.0358736139023148, 0.003165837697302747], gphase = [-0.8152414903269587, -2.782517169228026, 0.6712580617734987, -0.8740730716665552, -2.643571291758159, 0.47075081439120414, -0.9239140240255276, -2.5716056123269175, 0.27355110225263823, -0.9742703505846365  …  -0.5833398913691082, 0.13454898746871963, 1.5603672864839597, -0.24394569243170222, -2.938829994405641, -0.522251263172072, -0.02376437665701033, 1.5743124152515442, -0.4728503020642682, -2.9233844979917194]), (c = (params = [0.36746183235815133 0.067173405483948 … 0.6023681698716328 0.07503263231491032; 1.3611636643372957 1.1726998455330795 … 0.7832519939325823 0.7596371644008957; … ; -0.22362750845601753 -0.11205079397698516 … 1.122781240146348 1.5035678656988423; 0.96510489314312 -0.6820162896835006 … 0.7337305731674205 1.0811952937264973], hyperparams = 31.95335129999868), fg = 0.2940079986932449, σimg = 0.929384741374524, lgamp = [0.04611232661197692, -0.04705578600480994, -0.9809604706603325, -0.057727079221443475, 0.006396687260666782, -0.004306083407947614, -0.3495817663089692, 0.10836531223863463, -0.02012548159507765, 0.02021131090825903  …  -0.07728087352227744, 0.008911752720273279, -0.9807751246695559, 0.0001565908773751769, -0.05085238782238217, 0.04725317005026326, -0.054226400368089955, 0.00010418891743333727, -1.012269957548756, 0.00945355747334327], gphase = [-0.8133905870863778, -2.7261597680870038, 0.7167426251696539, -0.8712413736378788, -2.6433114925389445, 0.5088160403613923, -0.9257216600619972, -2.5410986123981387, 0.318555361968531, -0.9764983988581157  …  -0.5866809966145573, 0.033034942403883806, 1.459901430206287, -0.3049443640246806, -3.0659807841328344, -0.5223010544295057, -0.11830546461033313, 1.466324821324271, -0.5106349195869829, -3.028786339028994]), (c = (params = [0.21275783159633774 0.6407483090471261 … 0.9049995044084568 0.8903167335767875; -0.14239119202557374 -0.16903312881266236 … 0.857007589465344 0.4470440070379215; … ; 0.4568419060921419 -0.41903618464707726 … 1.116161741090328 -0.20658313962195482; -0.7512364169040825 0.27176656721175607 … 0.7741532662097029 -0.23736201867937703], hyperparams = 9.819108387074857), fg = 0.28061943090481367, σimg = 0.9788001289664419, lgamp = [0.022078908762314248, -0.01568837265066116, -0.8988630206597703, -0.002309044194348842, -0.0055522186704844755, 0.005824615490686925, -0.3246729338261919, 0.16232013261890899, -0.01902158546100462, 0.02535989599756591  …  -0.05588062051086281, 0.0029437999507065603, -0.9872584046831581, -0.015737311991191483, -0.0007708486900069511, 0.0002906315264964072, -0.10200614739592204, -0.01421144108718409, -1.0537534482974378, -0.005223490217451335], gphase = [-0.8168461988788529, -2.63187875541104, 0.6187139869178924, -0.8720863366440647, -2.5775061286306613, 0.42116561171095734, -0.930428252790174, -2.4872936816679028, 0.22830908646392814, -0.9771919931788209  …  -0.5866997500359716, 0.05897756933366417, 1.6603218625503748, -0.24174318023374414, -2.864750331879269, -0.5257359801682789, -0.060294149736137184, 1.6449651547359736, -0.45314018802165584, -2.8405375757097975]), (c = (params = [-0.704280822329594 -1.6671138774479142 … 0.10015460072629144 -0.20220416778132003; -0.4983711002030803 -0.6755693367421948 … 0.11575811263031155 -0.43964204563839443; … ; -0.3031220415128515 -0.4375518208507009 … 0.627396188008587 1.4945458424224924; 0.44638512601954566 -0.9953291756424755 … 1.00213072430607 1.4228243985819085], hyperparams = 23.156862629906264), fg = 0.3061040247230097, σimg = 0.9084095087601987, lgamp = [-0.03273162252322942, 0.034514268047885724, -0.9045333800584839, 0.0011341124678735208, -0.000727858075032728, 0.0004833581373707249, -0.3697802065323648, 0.12670868194186527, -0.018163197492546093, 0.02026671546947199  …  -0.059559225038693535, -0.002172789262284554, -1.0004444901409209, 0.0051132721272771665, -0.013157614293586816, 0.01422710925021428, -0.06666801372198476, 0.00689989141623629, -1.0519756336337995, 0.0045444944157427836], gphase = [-0.8192744871827039, -2.647216091538378, 0.5056695073277786, -0.8711946114843679, -2.5803823054996324, 0.2989443938142344, -0.9293993332223871, -2.503693171917364, 0.09945506810282392, -0.981340105131263  …  -0.5828980014937342, 0.14210701774341752, 1.6265439218673459, -0.24524037586036904, -2.8940931574807367, -0.5219295941717983, -0.016276449015789354, 1.6368657076176025, -0.4308079763233138, -2.8675999623785593]), (c = (params = [-0.8342933002417593 0.007210502091662345 … 0.45698641753330005 -0.3124659634446086; -0.883099481272986 -0.8236167826216712 … -0.623796503854381 -0.998843416591027; … ; 0.3635883012781486 -0.2946571633248756 … 0.7469195802915907 -0.009437212403762689; -0.753042999942317 -0.12858474525679908 … 0.2508195708712166 -0.6555903540396466], hyperparams = 6.043755349647896), fg = 0.2769322606031637, σimg = 0.9269500213561123, lgamp = [-0.027065579849175236, 0.030293393047732733, -0.9735636742742821, 0.033129461121535954, -0.0051034232157282, 0.00470263076121899, -0.4034500571229471, 0.1486656413826135, -0.027963815083049424, 0.02581525128855674  …  -0.06956913412671545, -0.019930237288314404, -0.9934878669026598, -0.009117591481974678, -0.009506395440056017, 0.01690594225596483, -0.057500993876432066, -0.005126733293260016, -1.039230640410123, 0.014713882700650598], gphase = [-0.8098933939431344, -2.7195808175434957, 0.5016456893663486, -0.8716817428915407, -2.6821571564283486, 0.28198904981435235, -0.9266585857164981, -2.5882660247894034, 0.10004242816637203, -0.9815624004322859  …  -0.5847591980396261, 0.11778373801733234, 1.5648360181075545, -0.21492766468137942, -2.9350426942848578, -0.5262825796931694, -0.030642212704945986, 1.5777266178750562, -0.4535014922809678, -2.925023872785013]), (c = (params = [-0.22778143958776506 0.35540783976571017 … -0.6558877211976675 0.32587542599173996; 0.07472152341249358 0.27622540025065506 … -0.9376688533668408 0.1290947697521241; … ; 0.02080132453881734 -0.845417413515031 … -0.20721082947293942 -0.3545366597712077; 0.11756278015984839 0.04189689741318059 … -0.26145032937358353 0.1395926621170575], hyperparams = 39.67996311152073), fg = 0.2768107252363007, σimg = 1.0074464302756836, lgamp = [-0.0009340440772038613, 0.0012812614445753714, -1.0528859506173087, -0.031461413485690176, -0.02693612904534547, 0.024807300314074462, -0.4407162830149158, 0.11586884442043276, -0.0156279327406478, 0.014216200397418808  …  -0.027443969088076517, 0.004469358394838764, -0.9625132273568304, -0.006569651028451951, -0.026605804190760415, 0.029396602818504525, -0.07579169142824177, 0.0027495185432368113, -1.0486082107739199, -0.009640196218376019], gphase = [-0.8102927431859639, -2.935074253819748, 0.6598821178455452, -0.8706253494662585, -2.7991342013956935, 0.47575494308939925, -0.9279139535659995, -2.6935714030117848, 0.2861231466896801, -0.978636456798987  …  -0.5831361110523731, 0.0014315021863040359, 1.4329057044845657, -0.3666923688104342, -3.0657630244388887, -0.5257413052253583, -0.13349417471053934, 1.4197679878349505, -0.5768474284773967, -3.0759149599680935]), (c = (params = [1.1423077912978 -0.3358877989453126 … -0.12159603283606446 -0.5301978526720256; 0.685498045583878 -0.50460198573896 … -0.2721442511552352 -0.022922301655938445; … ; -0.6734231312411872 -0.4393699556016573 … -1.4172052577885923 -0.4525166459487403; 0.03754652891334638 -0.40625688053582354 … -0.7307939451853811 -0.6294877956593984], hyperparams = 10.27021664294877), fg = 0.3137468171204755, σimg = 1.017152165305182, lgamp = [-0.01616495683564692, 0.01597794842969767, -1.081752465973352, -0.003496712304131073, -0.0205194705415981, 0.019347165126749617, -0.45078439811468374, 0.11640625470096058, -0.03478884912249393, 0.03709225698238195  …  -0.05353027683208388, 0.0021541074196171914, -0.9915944151214068, -0.004056938560137771, -0.05237957871971887, 0.05326547238635795, -0.04029118506392933, 0.006411582919922188, -1.0442216156360937, -0.008000584503451017], gphase = [-0.817266091602763, -2.940334030567806, 0.494620992035209, -0.8729075903831374, -2.826378550864346, 0.3055574275844306, -0.9300843179516157, -2.73724755393974, 0.11414514993640204, -0.9760898812113232  …  -0.5824171486334166, -0.004933774203091012, 1.4902507127481677, -0.38584982415116353, -3.030363745331529, -0.5295523189622571, -0.14895332898727892, 1.4834218222616362, -0.6177086807576784, -3.008920056938642])], NamedTuple{(:n_steps, :is_accept, :acceptance_rate, :log_density, :hamiltonian_energy, :hamiltonian_energy_error, :max_hamiltonian_energy_error, :tree_depth, :numerical_error, :step_size, :nom_step_size, :is_adapt), Tuple{Int64, Bool, Float64, Float64, Float64, Float64, Float64, Int64, Bool, Float64, Float64, Bool}}[(n_steps = 1023, is_accept = 1, acceptance_rate = 0.9927079457302517, log_density = 2936.7645144514445, hamiltonian_energy = -2385.7696126278825, hamiltonian_energy_error = 0.007928306648409489, max_hamiltonian_energy_error = 0.010703655803354195, tree_depth = 10, numerical_error = 0, step_size = 0.0001, nom_step_size = 0.0001, is_adapt = 1), (n_steps = 2, is_accept = 1, acceptance_rate = 1.0443148661842367e-36, log_density = 2936.7645144514445, hamiltonian_energy = -2282.832854659685, hamiltonian_energy_error = 0.0, max_hamiltonian_energy_error = 5039.920415122823, tree_depth = 1, numerical_error = 1, step_size = 0.0015547071097859802, nom_step_size = 0.0015547071097859802, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9442682403541073, log_density = 2517.980768566077, hamiltonian_energy = -2242.0945644208346, hamiltonian_energy_error = 0.08058731851588163, max_hamiltonian_energy_error = 0.11574905218776621, tree_depth = 10, numerical_error = 0, step_size = 0.0003024924868374592, nom_step_size = 0.0003024924868374592, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9462217055867638, log_density = 2271.46511695999, hamiltonian_energy = -1844.1783165157192, hamiltonian_energy_error = -0.02972853503456463, max_hamiltonian_energy_error = 0.184632359164425, tree_depth = 10, numerical_error = 0, step_size = 0.00043425969946870484, nom_step_size = 0.00043425969946870484, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.8644194766909357, log_density = 2060.8275142403813, hamiltonian_energy = -1636.9095258235275, hamiltonian_energy_error = 0.18660928839835833, max_hamiltonian_energy_error = 0.3773767960033183, tree_depth = 10, numerical_error = 0, step_size = 0.0007162540836695529, nom_step_size = 0.0007162540836695529, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.8413441925477443, log_density = 1887.2852289878915, hamiltonian_energy = -1405.0724441327125, hamiltonian_energy_error = -0.036078735156479524, max_hamiltonian_energy_error = 0.7552521506920584, tree_depth = 10, numerical_error = 0, step_size = 0.0009929215361010615, nom_step_size = 0.0009929215361010615, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.0005077388363158758, log_density = 1889.7595087373916, hamiltonian_energy = -1159.4851105538132, hamiltonian_energy_error = 0.6773073003782883, max_hamiltonian_energy_error = 256.4797369598333, tree_depth = 10, numerical_error = 0, step_size = 0.0013130964346827535, nom_step_size = 0.0013130964346827535, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9986654565353855, log_density = 1842.2056145849238, hamiltonian_energy = -1213.657299919872, hamiltonian_energy_error = 0.0023625062922292273, max_hamiltonian_energy_error = 0.007010090102539834, tree_depth = 10, numerical_error = 0, step_size = 0.00012796481185624653, nom_step_size = 0.00012796481185624653, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.995495654581789, log_density = 1849.8146779629838, hamiltonian_energy = -1201.6893108187644, hamiltonian_energy_error = -0.0033605213595819805, max_hamiltonian_energy_error = -0.039840237058797356, tree_depth = 10, numerical_error = 0, step_size = 0.00027407113605874394, nom_step_size = 0.00027407113605874394, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.8735129349668224, log_density = 1646.8011149143038, hamiltonian_energy = -1168.7104992567424, hamiltonian_energy_error = 0.16094071125962728, max_hamiltonian_energy_error = 0.40359894958601217, tree_depth = 10, numerical_error = 0, step_size = 0.0005913320057177703, nom_step_size = 0.0005913320057177703, is_adapt = 1)  …  (n_steps = 1023, is_accept = 1, acceptance_rate = 0.8072608317400083, log_density = 1169.2603779932779, hamiltonian_energy = -493.71717223088376, hamiltonian_energy_error = 0.14235355918572168, max_hamiltonian_energy_error = 0.6804148696897983, tree_depth = 10, numerical_error = 0, step_size = 0.004730444126673395, nom_step_size = 0.004730444126673395, is_adapt = 0), (n_steps = 511, is_accept = 1, acceptance_rate = 0.9370451539897652, log_density = 1130.4354465872184, hamiltonian_energy = -470.76936744363184, hamiltonian_energy_error = 0.09803598950657033, max_hamiltonian_energy_error = 0.2784770948355799, tree_depth = 9, numerical_error = 0, step_size = 0.004730444126673395, nom_step_size = 0.004730444126673395, is_adapt = 0), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.8968105046582765, log_density = 1119.3071081267083, hamiltonian_energy = -461.39423345247985, hamiltonian_energy_error = -0.3609665373445523, max_hamiltonian_energy_error = 0.6312017135678616, tree_depth = 10, numerical_error = 0, step_size = 0.004730444126673395, nom_step_size = 0.004730444126673395, is_adapt = 0), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.6500482394886792, log_density = 1138.2320602757065, hamiltonian_energy = -446.3091057903091, hamiltonian_energy_error = 0.2668942276548023, max_hamiltonian_energy_error = 1.3556924371889636, tree_depth = 10, numerical_error = 0, step_size = 0.004730444126673395, nom_step_size = 0.004730444126673395, is_adapt = 0), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9623765359964752, log_density = 1138.5588715940655, hamiltonian_energy = -442.3434376659818, hamiltonian_energy_error = -0.003803624175816367, max_hamiltonian_energy_error = -0.6380739355206515, tree_depth = 10, numerical_error = 0, step_size = 0.004730444126673395, nom_step_size = 0.004730444126673395, is_adapt = 0), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9849251207219568, log_density = 1159.1964927630013, hamiltonian_energy = -465.44533102229957, hamiltonian_energy_error = -0.1532954414062715, max_hamiltonian_energy_error = -0.7028802645832002, tree_depth = 10, numerical_error = 0, step_size = 0.004730444126673395, nom_step_size = 0.004730444126673395, is_adapt = 0), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.8063614156425226, log_density = 1144.3844449378917, hamiltonian_energy = -476.34671660432457, hamiltonian_energy_error = -0.05960899601984693, max_hamiltonian_energy_error = 0.911753970106929, tree_depth = 10, numerical_error = 0, step_size = 0.004730444126673395, nom_step_size = 0.004730444126673395, is_adapt = 0), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.882288152227988, log_density = 1141.6510178938888, hamiltonian_energy = -454.8528412719012, hamiltonian_energy_error = 0.12465676900376366, max_hamiltonian_energy_error = 0.668090328200833, tree_depth = 9, numerical_error = 0, step_size = 0.004730444126673395, nom_step_size = 0.004730444126673395, is_adapt = 0), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.8242349724107931, log_density = 1121.6293259670758, hamiltonian_energy = -443.20627339541056, hamiltonian_energy_error = 0.3226021256317608, max_hamiltonian_energy_error = 1.1044954766394994, tree_depth = 10, numerical_error = 0, step_size = 0.004730444126673395, nom_step_size = 0.004730444126673395, is_adapt = 0), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9713618718513379, log_density = 1136.489539351714, hamiltonian_energy = -418.40587737723695, hamiltonian_energy_error = -0.046441054694696504, max_hamiltonian_energy_error = -1.314191565907322, tree_depth = 10, numerical_error = 0, step_size = 0.004730444126673395, nom_step_size = 0.004730444126673395, is_adapt = 0)])
Note

The above sampler will store the samples in memory, i.e. RAM. For large models this can lead to out-of-memory issues. To fix that you can include the keyword argument saveto = DiskStore() which periodically saves the samples to disk limiting memory useage. You can load the chain using load_table(diskout) where diskout is the object returned from sample. For more information please see ComradeAHMC.

Now we prune the adaptation phase

chain = chain[501:end]
Table with 5 columns and 200 rows:
      c                     fg        σimg      lgamp                 ⋯
    ┌──────────────────────────────────────────────────────────────────
 1  │ (params = [-0.62974…  0.274379  1.03642   [0.00572853, -0.006…  ⋯
 2  │ (params = [0.090859…  0.287552  1.03405   [0.00407881, -0.007…  ⋯
 3  │ (params = [-0.55023…  0.281611  0.973517  [-0.00659724, 0.003…  ⋯
 4  │ (params = [-0.29607…  0.328127  1.03077   [-9.27783e-6, -0.00…  ⋯
 5  │ (params = [-0.86711…  0.307151  0.96579   [-0.04786, 0.045038…  ⋯
 6  │ (params = [-0.66657…  0.299196  0.915323  [-0.0304487, 0.0262…  ⋯
 7  │ (params = [-1.10293…  0.30599   0.93328   [-0.0136521, 0.0170…  ⋯
 8  │ (params = [0.224554…  0.355963  0.972824  [-0.0160305, 0.0141…  ⋯
 9  │ (params = [-0.43721…  0.27889   0.883478  [-0.0186111, 0.0215…  ⋯
 10 │ (params = [-0.10894…  0.279788  0.846361  [-0.012547, 0.01085…  ⋯
 11 │ (params = [-0.21589…  0.291033  0.880974  [-0.00830623, 0.007…  ⋯
 12 │ (params = [0.242064…  0.308047  0.942483  [-0.0535764, 0.0514…  ⋯
 13 │ (params = [0.029459…  0.310245  0.851606  [-0.00901956, 0.008…  ⋯
 14 │ (params = [-1.41495…  0.33107   0.899375  [-0.0100774, 0.0158…  ⋯
 15 │ (params = [0.426702…  0.313856  0.931994  [-0.0321854, 0.0339…  ⋯
 16 │ (params = [-1.71619…  0.316186  0.984115  [-0.0317408, 0.0274…  ⋯
 17 │ (params = [-1.44036…  0.311323  0.957748  [-0.0309214, 0.0301…  ⋯
 ⋮  │          ⋮               ⋮         ⋮               ⋮            ⋱
Warning

This should be run for likely an order of magnitude more steps to properly estimate expectations of the posterior

Now that we have our posterior, we can put error bars on all of our plots above. Let's start by finding the mean and standard deviation of the gain phases

gphase  = hcat(chain.gphase...)
mgphase = mean(gphase, dims=2)
sgphase = std(gphase, dims=2)
104×1 Matrix{Float64}:
 0.003557682417846069
 2.1031427549678012
 0.44684804960222385
 0.003060687163304964
 1.5073793300889649
 0.4510136990744012
 0.003069702269742162
 0.7423338716868217
 0.45390721972401904
 0.003094223438640786
 ⋮
 0.28156599561706575
 0.400445525604335
 0.2239061886675918
 2.7906478929332486
 0.003922085060623269
 0.27703278746727855
 0.3957879490839807
 0.21879937333661115
 2.8018964719353807

and now the gain amplitudes

gamp  = exp.(hcat(chain.lgamp...))
mgamp = mean(gamp, dims=2)
sgamp = std(gamp, dims=2)
129×1 Matrix{Float64}:
 0.02116149483471646
 0.021340142650787236
 0.027099673801421212
 0.03677923989663481
 0.01717261463750595
 0.017207672279772736
 0.03805269681212777
 0.03926505580049912
 0.0170764805946439
 0.017853763083506554
 ⋮
 0.009049293314805612
 0.015392921806148938
 0.009020437572396465
 0.015805455238790064
 0.016251252634688327
 0.02457868690735038
 0.010098542010158685
 0.015145597791362628
 0.009669335170843632

Now we can use the measurements package to automatically plot everything with error bars. First we create a caltable the same way but making sure all of our variables have errors attached to them.

using Measurements
gmeas_am = measurement.(mgamp, sgamp)
ctable_am = caltable(gcache, vec(gmeas_am)) # caltable expects gmeas_am to be a Vector
gmeas_ph = measurement.(mgphase, sgphase)
ctable_ph = caltable(gcachep, vec(gmeas_ph))
───────────┬────────────────────────────────────────────────────────────────────
      time │      AA            AP         AZ          JC          LM          ⋯
───────────┼────────────────────────────────────────────────────────────────────
 0.917±0.0 │ 1.0±0.0  -0.815±0.004    missing     missing    -1.9±2.1   0.58±0 ⋯
 1.217±0.0 │ 1.0±0.0  -0.871±0.003    missing     missing    -2.3±1.5   0.38±0 ⋯
 1.517±0.0 │ 1.0±0.0  -0.928±0.003    missing     missing  -2.59±0.74   0.18±0 ⋯
 1.817±0.0 │ 1.0±0.0  -0.977±0.003    missing     missing  -2.56±0.62  0.009±0 ⋯
 2.117±0.0 │ 1.0±0.0  -1.003±0.004    missing     missing  -2.55±0.25  -0.16±0 ⋯
  2.45±0.0 │ missing       1.0±0.0    missing     missing    0.66±2.9   0.65±0 ⋯
  2.75±0.0 │ 1.0±0.0  -1.046±0.007    missing     missing  -2.48±0.25  -0.54±0 ⋯
  3.05±0.0 │ 1.0±0.0  -1.132±0.003    missing     missing  -2.44±0.25  -0.72±0 ⋯
  3.35±0.0 │ 1.0±0.0  -1.152±0.003    missing     missing  -2.44±0.25  -0.89±0 ⋯
 3.683±0.0 │ 1.0±0.0  -1.163±0.005   0.77±2.8     missing   -2.5±0.25  -1.08±0 ⋯
 3.983±0.0 │ 1.0±0.0  -1.167±0.004    1.7±2.2     missing  -2.54±0.25  -1.25±0 ⋯
 4.283±0.0 │ 1.0±0.0   -1.16±0.004    2.3±1.4     missing  -2.58±0.47  -1.38±0 ⋯
 4.583±0.0 │ 1.0±0.0  -1.143±0.003  2.47±0.63  -0.22±0.37    -2.5±1.2  -1.56±0 ⋯
 4.917±0.0 │ 1.0±0.0  -1.073±0.006   2.33±0.3  0.038±0.38    -1.8±2.2  -1.74±0 ⋯
 5.183±0.0 │ 1.0±0.0  -1.078±0.006   2.17±0.3   0.22±0.39   -0.84±2.8  -1.88±0 ⋯
  5.45±0.0 │ 1.0±0.0   -1.119±0.01   2.03±0.3    0.41±0.4     1.1±2.7  -2.01±0 ⋯
     ⋮     │    ⋮          ⋮            ⋮          ⋮           ⋮           ⋮   ⋱
───────────┴────────────────────────────────────────────────────────────────────
                                                    2 columns and 9 rows omitted

Now let's plot the phase curves

plot(ctable_ph, layout=(3,3), size=(600,500))
Example block output

and now the amplitude curves

plot(ctable_am, layout=(3,3), size=(600,500))
Example block output

Finally let's construct some representative image reconstructions.

samples = skymodel.(Ref(post), chain[begin:2:end])
imgs = intensitymap.(samples, fovx, fovy, 128,  128)

mimg = mean(imgs)
simg = std(imgs)
fig = CM.Figure(;resolution=(400, 400))
CM.image(fig[1,1], mimg,
                   axis=(xreversed=true, aspect=1, title="Mean Image"),
                   colormap=:afmhot)
CM.image(fig[1,2], simg./(max.(mimg, 1e-5)),
                   axis=(xreversed=true, aspect=1, title="1/SNR",),
                   colormap=:afmhot)
CM.image(fig[2,1], imgs[1],
                   axis=(xreversed=true, aspect=1,title="Draw 1"),
                   colormap=:afmhot)
CM.image(fig[2,2], imgs[end],
                   axis=(xreversed=true, aspect=1,title="Draw 2"),
                   colormap=:afmhot)
fig

Now let's check the residuals

p = plot();
for s in sample(chain, 10)
    residual!(p, vlbimodel(post, s), dvis)
end
p
Example block output

And viola, you have just finished making a preliminary image and instrument model reconstruction. In reality, you should run the sample step for many more MCMC steps to get a reliable estimate for the reconstructed image and instrument model parameters.

This page was generated using Literate.jl.